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I.  INTRODUCTION 


A.  Scope 

This  report  summarizes  the  results  of  the  initial  phase  of  a  comprehen¬ 
sive  simulation  study  of  alternative  signal  processing  algorithms  for  data 
adaptive  superresolution  direction  finding  and  spatial  nulling  to  support  sig¬ 
nal  copy  in  the  presence  of  strong  cochannel  Interference.  The  need  for  such 
a  study  arises  because,  although  most  of  the  techniques  evaluated  have  been 
documented  in  the  literature,  no  systematic  comparison  has  heretofore  been 

undertaken. 

The  general  approach  of  the  current  study  is  to  simulate  a  sequence  of 
increasingly  more  general,  i.e,  realistic,  signaling  environments,  and  to 
expose  each  of  the  more  promising  algorithms  to  all  of  the  "standardized" 
experiments,  in  turn.  For  the  initial  phase  of  the  inquiry,  we  have  selected 
an  Ideal  environment  characterized  by  a  uniform,  linear  array  of  identical 
isotropic  elements  and  perfect  receivers.  In  addition,  partial  results  are 
obtained  for  the  case  when  the  array  steering  vectors  are  in  error  by  a  small 
amount  which  might  be  caused  by  residual  calibration  errors,  unmodeled  multi- 
path  distortions,  near-field  emitters,  etc. 

B.  Focus  of  the  Initial  Inquiry 

Many  of  the  techniques  proposed  for  superresolution  array  processing  have 
their  origin  in  spectral  estimation  for  time  series.  Since  the  sampling  of  a 
function  in  time  is  analogous  to  sampling  a  function  in  space,  it  is  natural 
to  make  this  association;  estimating  the  frequency  of  sinusoids  in  noise  can 
be  seen  to  be  equivalent  to  estimating  the  directions  of  planewaves  in  noise. 

Although  the  superresolution  problem  involves  finding  plane  waves  in 
noise,  most  spectral  estimation  techniques  make  no  use  of  any  information 
about  the  underlying  process.  Indeed,  such  methods  are  only  heuristically 
motivated,  since  the  estimation  of  a  completely  unknown  function  based  upon  a 
finite  number  of  samples  is  at  best  an  underdefined  process.  Nevertheless, 
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since  claims  of  success  have  been  made  for  such  techniques,  we  thought  it  best 
to  start  our  investigation  by  reviewing  these  "classical”  techniques. 
5.epresentative  of  such  techniques  are  the  following: 

(1)  Adapted  Angular  Response  (AAR)  [13] 

(2)  Maximum  Entropy  Method  (MEM)  [9] 

(3)  Maximum  Likelihood  Method  (MLM)  [12] 

(4)  Thermal  Noise  Algorithm  (TNA)  [15] 

By  contrast  with  classical  spectral  methods,  a  technique  which  uses  (re¬ 
quires)  the  fact  that  the  number  of  plane  waves  is  finite  is  the  MUSIC 
algorithm  [16],  MUSIC,  which  denotes  Multiple  Signal  £lassificatlon,  is  an 
extension  of  the  method  of  Pisarenko  [18].  MUSIC  is  but  one  member  of  a  class 
of  methods  based  upon  the  decomposition  of  covariance  data  into  eigenvectors 
and  eigenvalues.  Such  techniques,  known  as  Singular  Value  Decompositions 
(SVD’s),  will  be  more  completely  reviewed  in  a  subsequent  phase  of  the  study. 
In  order  to  place  SVD  techniques  relative  to  the  classical  methods,  however, 
results  for  MUSIC  are  included  in  this  report. 

All  of  the  techniques  reviewed  have  application  to  arbitrary  array  geom¬ 
etries.  When  the  array  happens  to  be  linear,  however,  it  is  possible  to  take 
advantage  of  the  structure  of  the  array  to  Improve  the  resolving  power  of 
these  techniques  [2].  Each  of  the  techniques  involves  the  calculation  of  a 
quadratic  form  which,  in  the  case  of  the  linear  array,  becomes  a  complex, 
trigonometric  polynomial  which  is  being  evaluated  on  the  unit  circle  in  the 
complex  plane.  The  behavior  of  such  functions  is  dominated  by  the  polynomial 
roots  which  lie  on  or  close  to  the  unit  circle.  Thus,  a  rooting  variant  of 
each  technique  is  possible  which  computes  the  arguments  of  the  roots  closest 
to  the  unit  circle.  The  rooting  variant  of  TNA  is  identical  to  that  of  AAR, 
so  that  only  four  rooting  algorithms  have  been  presented. 
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The  adaptive  listening  or  copy  functions  for  closely-spaced  signals 
requires  the  use  of  main-beam  nulling.  Since  this  is  difficult  to  accomplish 
without  accidentally  nulling  the  desired  signal  [23],  however,  much  attention 
must  he  placed  upon  the  sensitivity  introduced  by  only  approximately  knowing 
the  direction  of  the  desired  signal.  Two  covariance  modeling  techniques  for 
steering  nulls  in  the  direction  of  Interferers  are  reviewed  in  this  report. 

The  sensitivity  to  errors  in  the  knowledge  of  the  array  response  to  plane 
waves  will  be  the  subject  of  a  detailed  study  in  a  subsequent  report.  Only  a 
brief  introduction  to  this  subject  is  included  here  to  determine  the  conse¬ 
quences  of  such  errors  on  the  ability  to  place  nulls  for  copy  by  the  modeling 
techniques  described  above.  A  simple  direction-independent  gain  and  phase 
error,  independent  from  receiver  to  receiver,  has  been  utilized  for  this 
initial  investigation. 

C.  Overview  of  the  Report 

The  following  section  (II)  describes  the  performance  bounds  for  direction 
finding  and  copy.  The  bounds  for  estimation  errors  are  the  Cramer-Rao  type 
for  the  assumed  case  of  Gaussian  signals  in  Gaussian  noise;  the  limiting  per¬ 
formance  for  ideal  arrays  and  infinite  data  are  given  for  copy  performance 
bounds • 

An  expanded  description  of  the  algorithms  considered  during  this  study  is 
given  in  Section  III.  The  emphasis  in  this  section  is  a  comparative  introduc¬ 
tion  to  the  techniques  for  direction  finding,  power  estimation  and  adaptive 
copy  weighting. 

An  overview  of  the  simulation  approach  used  for  this  study  is  given  in 
Section  IV  and  the  results  are  summarized  in  Section  V.  Various  appendices 
elaborating  on  some  topics  and  a  complete  collection  of  the  Monte  Carlo  simu¬ 
lation  outputs  are  also  included. 
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D.  Major  Findings  ' 

The  principal  conclusions  based  upon  the  initial  experiments 
follows: 

1.  Although  some  of  the  spectral  algorithms  tested  to  date  are  more 
sensitive  than  MUSIC,  the  angle  estimates  provided  by  these 
algorithms  are  generally  poor. 

2.  For  the  linear  array  problem,  a  root  variant  of  MUSIC  exists 
which  is  considerably  more  sensitive  than  the  spectral  version. 

3.  The  direct  solution  for  the  power  present  in  multiple  directions 
can  fail  because  of  mutual  dependencies  among  direction  vectors 
whenever  the  array  is  irregular.  This  difficulty  can  be  overcome 
by  formulating  the  power  estimation  as  a  least-squares  problem. 

4.  Modeling  of  interference  vectors  provides  an  effective  means  of 
avoiding  signal  cancellation  as  a  result  of  direction— of —arrival 
errors  for  the  desired  signal.  Unfortunately,  the  resulting  per¬ 
formance  is  very  sensitive  to  array  calibration  errors. 

The  implications  of  these  observations  are  that  Singular-Value  Decomposi¬ 
tion  is  a  desirable  prerequisite  for  superresolution  of  plane  waves  and 
the  exploitation  of  special  array  structures  leads  to  increased  resolution 
sensitivity,  but  that  sensitivity  to  array  errors  is  the  most  serious  obstacle 
to  successful  implementation.  For  these  reasons,  the  emphasis  for  the  next 
phase  of  the  study  will  be  in  the  following  areas: 

1.  Singular-Value  Decomposition  techniques 

2.  Maximal  exploitation  of  linear  array  structures 

3.  Robust  algorithms  to  reduce  array  error  sensitivities 

In  addition  to  these  areas  of  emphasis,  Monte  Carlo  experiments  will  be 
extended  to  explore  the  effects  of  larger  numbers  of  signals,  decreased  SIR, 
and  various  forms  of  array  errors . 
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II.  THEORETICAL  PERFORMANCE  BOUNDS 


A.  Direction-Finding 

A  major  portion  of  this  report  is  concerned  with  direction-finding  (DF) 
algorithms  and  their  expected  performance,  as  determined  by  extensive  Monte 
Carlo  simulations.  Although  these  results  can  be  studied  on  a  stand-alone 
basis,  suitable  theoretical  benchmarks  serve  to  put  the  results  in  a  better 

perspective* 

Perhaps  the  most  desirable  benchmark  would  be  the  performance  achieved  by 
the  quintessential  ”optimum"  processor.  Given  such  a  benchmark,  one  could 
reasonably  expect  to  make  sound-  judgements  as  to  whether  or  not  a  particular 
algorithm  was  “good  enough".  Unfortunately,  optimum  processors  can  be  prohib¬ 
itively  expensive  to  simulate.  Moreover,  practical  applications  are  often 
sufficiently  complicated  that  a  single  criterion  for  "optimality"  is  virtually 
impossible  to  define* 

As  a  relevant  case  in  point,  consider  the  problem  of  direction-finding  in 
a  multiple  emitter  environment.  This  problem  is  characterized  by  an  unknown 
number  of  signals  arriving  from  unknown  directions.  Thus,  a  good  DF  algorithm 
must  determine  the  number  of  signals  present  as  well  as  provide  accurate  di¬ 
rection  estimates.  In  some  instances,  estimates  of  the  signal  power  levels 
are  also  required.  Finally,  in  order  to  be  useful  in  unfriendly  environments, 
all  of  these  requirements  must  be  met  without  detailed  knowledge  of  the  signal 
wavef  oirms* 

Taken  in  its  entirety,  the  multiple  emitter  DF  problem  is  too  complex  to 
admit  a  comprehensive  theoretical  analysis.  One  possible  simplification  is  to 
specify  the  number  of  emitters.  In  this  case,  one  can  obtain  theoretical 
bounds  on  the  accuracy  of  the  (DF)  estimates  by  computing  the  relevant  Fisher 
information  matrix.  Inverting  this  matrix  yields  the  Cramer-Rao  bound  on  the 
variance  of  any  unbiased  estitnate* 

For  the  purposes  of  direction-finding,  the  assumption  is  made  that  the 
signal  sources  are  stationary  (complex)  Gaussian  random  processes.  The  re- 
eived  signals  are  sampled  at  a  rate  less  than  the  receiver  bandwidth.  Under 


5 


the  latter  assumption,  the  signals  obtained  at  two  different  instants  of  time 
are  statistically  uncorrelated.  This  model  is  extremely  convenient  for  gener¬ 
ating  simulation  data  and  leads  to  performance  bounds  that  depend  only  on  the 
signal  directions  and  powers# 

The  Cramer-Rao  (C-R)  bound  for  locating  Gaussian  emitters  is  derived  in 
Appendix  A,  Here,  we  merely  state  the  underlying  model  and  present  the  final 
form  of  the  result.  To  check  the  tightness  of  the  bound,  F.  White  [1]  has 
investigated  the  performance  of  two  "optimum"  DF  processors.  His  results 
indicate  that  the  C-R  bound  is  achievable  over  an  interesting  and  broad  range 
of  parameter  values# 

The  samples  obtained  from  the  array  elements  at  any  given  instant  of 
time,  called  a  snapshot,  may  be  modelled  as  a  vector 

r  =  s  +  n  , 

where  n  denotes  the  contribution  from  thermal  noise  and  other  sources  of  error 
in  the  receiver(s)  and  s  represents  the  received  slgnal(s).  When  only  one 
emitter  is  present,  the  signal  vector  may  be  written  as 

s  =  v(a)  p  , 

where  p  is  the  complex  amplitude  of  the  signal  that  would  be  observed  at  the 
phase  center  of  the  array,  and  v(a)  is  a  vector  constructed  from  the  complex 
oltage  gains  of  the  individual  array  elements.  The  vector  v(a)  is  often 
referred  to  as  a  "direction"  or  "steering"  vector  since  it  depends  only  upon 
the  direction  of  arrival  of  the  signal.  For  linear  arrays,  the  most 
convenient  measure  of  direction  is  the  cosine  of  the  angle  between  the  array 
axis  and  the  llne-of-sight  to  the  signal.  The  direction  cosine  is  denoted  by 

a# 

When  more  than  one  emitter  is  present,  the  principle  of  superposition 
allows  us  to  write  the  received  signal  as 


6 


s  =  Vp 
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where  the  jth  element  of  the  vector  p  is  the  complex  amplitude  of  the  signal 
from  the  jth  emitter  as  "seen”  at  the  array  phase  center.  Note  that  the  jth 
column  of  the  matrix  V  is  the  direction  vector  associated  with  the  jth 
einitter,  i.e., 

V  =  [v(aj^)  j  •  •  •  j 

where  J  denotes  the  actual  number  of  emitters. 

Under  our  statistical  assumptions,  the  ”signal-in-space‘  vector  p  is 
completely  described  by  its  mean  value,  assumed  to  be  zero,  and  its  covariance 


P  =  E{pp^}  , 

where  E{x}  generally  denotes  the  expected  value  of  x,  and  a  superscript  H 
indicates  the  conjugate  (Hermltlan)  transpose  operation.  Naturally,  the 
Gaussian  model  is  extended  to  include  the  noise  vector  n,  and  we  represent  the 
noise  covariance  matrix  as 

N  =  E  {nn^  } 

Assuming  the  noise  statistics  are  known,  one  could  always  normalize  (or 
transform)  the  data  in  such  a  way  that  the  components  of  the  noise  vector  n 
are  identically  distributed  and  statistically  independent  (uncorrelated). 
Unless  otherwise  explicitly  stated,  we  will  proceed  under  the  assumption  that 
the  noise  covariance  is  the  identity  matrix  I. 

Consider  the  problem  of  estimating  the  directions  of  arrival  when  the 
signal-in-space  covariance  P  is  known.  The  Fisher  information  matrix  for  this 
problem  can  be  stated  in  a  reasonably  compact  form  by  first  introducing  the 
gramian  matrix 
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w  =  v®v 


and  solving 

P-0  =  PWO 

for  0.  It  can  be  shown  that  this  equation  always  has  a  unique  (Hermician) 
solution.  We  next  introduce  the  "derivative”  of  V 

V  =  [vCoj)  I  .  .  ♦  I  v(aj)]  , 

where 


V  =  dv/do 

is  the  usual  derivative  of  v  with  respect  to  o.  It  is  also  convenient  to 
define 

w  i  V®V 

Under  the  conditions  stated  above,  the  Fisher  matrix  for  the  unknown 
directions  a  is  given  by 

F  =  2K  Re  {  (P-0)’^  □  (v“v  -  VpQW)  +  (OW)^O  (QW)  } 
aa 

where  K  is  the  number  of  available  snapshots  (observations).  Re{x}  denotes 
the  real  part  of  x,  and  a  superscript  "T"  refers  to  the  usual  transpose 
operation.  The  element  by  element  (Hadamard)  product  of  two  matrices  A  and  B 
with  the  same  dimensions  is  written  as  AQB* 
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The  nain  diagonal  of  the  inverse  of  the  Fisher  matrix  provides  a  lower 
bound  on  the  accuracy  (i.e.,  variance)  of  any  unbiased  DF  estimate.  However, 
when  the  signal-in-space  covariance  P  is  also  unknown,  the  resulting  bounds 
are  not  as  tight  as  they  might  be. 

If  the  signal-in-space  covariance  were  completely  unspecified,  generating 
the  required  Fisher  matrix  would  become  extremely  awkward.  Fortunately,  many 
applications  of  interest  are  adequately  modelled  by  assuming  uncorrelated 
emitters.  In  this  important  special  case,  P  is  a  diagonal  matrix,  i.e., 

Plj  =  0  ;  i  9^  j 

Thus,  consider  the  vector  S  obtained  by  taking  the  logarithm  of  the  emitter 
"powers"  (i.e.,  the  main  diagonal  elements  of  P).  When  the  directions  of  ar¬ 
rival  are  known,  the  normalized  Fisher  matrix  for  estimating  the  emitter  pow- 
ers  can  be  written  as 

y  =  K  (OW)  □  (QW)"^  . 

PP 

When  the  directions  and  powers  are  both  unknown,  the  "reduced  information 
matrix  for  the  directions  is  generally  of  the  form 

=  F  -  F^  f”^  F  . 

aa  aoL  Sot  SB  Sot 

The  amount  of  information  lost  depends  upon  the  coupling  matrix,  which  is 
given  by 

F„  =  2K  Re  {  (OW)  □  (QW)’’’} 

6a 

for  the  problem  considered  here. 
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Inverting  the  reduced  information  matrix  yields  the  C-R  bound  for 
(unbiased)  DF  estimates  ohen  the  emitter  pooers  are  also  unknomn.  While  the 
exact  expressions  given  above  appear  to  be  quite  formidable,  the  asymptotic 
form  of  the  DF  bound  is  really  quite  simple.  As  the  signal-to-nolse  ratios  of 
all  of  the  emitters  become  arbitrarily  large,  me  may  replace  0  with  the  in¬ 
verse  of  W  (provided  it  exists,  of  course).  Eventually,  P  completely  domi- 
nates  0  and  the  asymptotic  approximation 

_>  F  — >  2K  ?QA 
aa  oto 

emerges,  where  the  array  factor 
A  = 

depends  only  on  the  directions  of  arrival.  Since  P  is  a  diagonal  matrix,  in¬ 
verting  FtlA  is  trivial  and  the  asymptotic  accuracy  of  any  unbiased  DF  esti- 

mate  is  bounded  by 


Var  {a^}  — >  [2K  P.  .A.^C ....  a^)] 

where  Var  {x}  denotes  the  variance  of  a  random  variable  x.  Surprisingly,  the 
asymptotic  error  predicted  by  the  DF  bound  is  independent  of  the  relative 
strength  of  the  emitters!  However,  as  expected,  the  asymptotic  DF  variance  is 
inversely  proportional  to  the  product  of  the  signal-to-noise  ratio  and  the 
number  of  snapshots • 

Unfortunately,  the  array  factor  is  generally  a  very  complicated  function 
of  the  directions  of  arrival  and  the  array  geometry.  However,  certain  simpli¬ 
fications  are  possible  under  the  assumption  of  identical  elements  in  a  linear 
array.  In  particular,  the  phase  reference  point  may  be  chosen  so  that 


v^v  =  0 

Consequently,  the  array  factor  for  a  single  emitter  is  given  by 
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and  the  array  factor  for  two  emitters  may  be  written  as 
02)  =  =  h(a2  -  ct^^)  ;  j  = 

where  n(a)  is  a  symmetric  function  which  approaches  unity  for  sufficiently 
large  a.  Thus,  n  eay  be  leterpreted  as  an  efficiency  factor  that  only  depends 
on  the  separation  of  the  two  emitters  [2]. 

Be  Adaptive  Listening 

A  phased  array  receiver  is  assumed  to  be  used  to  monitor,  or  listen,  to 
a  desired  emitter  by  discriminating  against  the  undesired  emitters  only  on  the 
basis  of  differences  in  angular  directions  of  arrival.  An  adaptive  array 
differs  from  a  more  conventional  phased  array  in  that  the  complex  weights 
associated  with  the  antenna  elements  are  not  determined  by  the  designer  a 
priori.  Instead,  these  weights  are  optimized,  according  to  some  criterion,  on 
the  basis  of  measurements  which  are  made  on  the  signal  environment. 

Let  xk  denote  vector  sample  of  complex  array  data  collected  during  the 
kth  snapshot,  where  k=l,...,K.  As  in  the  above  discussion,  the  data  are 

representable  as 

^  =  Ps^(«s^  Vk  ^  \ 

where  it  is  assumed  that  the  desired  emitter  arrives  from  direction  Og  with 
complex  amplitude  Pg  and  the  jth  component  of  the  vector  is  the  complex 
amplitude  of  the  jth  interference  wave  source  and 

Vi  =  [v(aj)|  .  .  . 

is  the  matrix  of  the  interferer  direction  vectors.  The  objective  of  the 
adaptive  listening  array  is  to  enhance  the  gain  of  the  array  toward  the 
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desired  signal  while  simultaneously  nulling  the  signal  energy  due  to  the 
directional  interference.  This  is  accomplished  by  computing  a  set  of 
complex-valued  element  weights  which,  when  applied  to  the  element  outputs  and 
linearly  combined,  results  in  a  complex  scalar  array  output  of  the  form 

H 

If  the  signals  and  noise  are  Gaussian,  the  weight  vector  which  maximizes  the 
probability  of  detection  of  the  desired  signal  is  given  by 

w  iC^  v(a  ) 
o  N  s 

where  %  is  the  spatial  covariance  of  the  interference  plus  noise,  viz. 

"n  =  '('''A  ‘Vk  * 


where,  as  above,  we  have  taken  the  noise  to  be  isotropic  with  unit  variance. 

Moreover,  this  same  weight  vector  maximizes  the  output  signal  power  to 
average  interference-plus-noise  power  ratio  (SIR),  given  by 


SIR  =  Pg|w\(a)|  /w”  W  , 

Where  is  the  power  of  the  desired  signal,  thus,  the  maximum  output  SIR  is 

s 

given  by 


SIR  =  v’^Ca  )  R^^ 
max  s  s  N 


vCa^) 
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In  the  present  application,  it  is  important  to  separate  the  desired 
emitter  from  the  interference  and  noise.  In  this  .case,  we  start  by  observing 
the  covariance  with  all  signals  present 


R 


Of  course,  in  the  limit  for  large  sample  size,  R  converges  to  R,  where 


v(o  )  v®(a  )■ 
s  s 


so  that 

R'^vCOg)  =  [1/(1 


where  T 

o 


SIR 


max 


Thus,  if  v(a  )  is  correct,  w  =  r"^  v(a  )  will  converge  to  the  maximum  SIR 
>  g  s 

weight  vector  in  the  limit  for  large  R.  Unfortunately,  as  Miller  [49]  and 
Boroson  [50]  have  noted,  the  convergence  rate  will  be  dependent  upon  Tq. 
Specifically,  if  we  define  a  generalized  SIR  as 


GIR(Vg,P,  v)  =  Pg 


V 

s 


V 


9 


then  the  convergence  of  the  sampled  data  weight  vector  depends  upon 


•^1 


A 


GIR(Vg,  R,  Vg) 
T 

o 


p/[l  +  (1  -  p)] 


where  p  is  distributed  according  to  the  Beta  density  [49],  [50] 
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f(p)  =  [1  /  B(L  -  1,  K  +  2  -  1)]  (1  -  P)^  ^  > 

0  <  p  <  1 

with  B(M,  N)  =  (M-1)! (N-1) ! 

Therefore , 

P[pj  _<  1  -  5]  =  P{p  1  (1  -  5)  (1  +  ^ 

If  the  desired  signal  power  in  the  covariance  matrix  were  zero,  the  above 
expression  predicts  the  familiar  result  of  Reed,  et  al.  [22]  that  K  2  2L  -  3 
will  result  in  Et pi]  2  ^''2.  On  the  other  hand,  with  Tq  =  20  dB,  as  many 
as  50  times  as  many  samples  would  be  needed  to  achieve  the  same  result  [50, 

Fig.  3]. 

Throughout  the  preceding  discussion,  the  steering  constraint  for  the 
listening  weight  vector  was  taken  to  be  correct.  Actually,  this  constraint 
must  be  obtained  from  the  direction-finding  process  described  above  in  Section 
A.  Thus,  we  have  only  an  estimate  Og  of  the  signal  direction.  Since  we 
need  a  direction  vector  to  constrain  the  listening  weights,  however,  we  may 
choose  to  use  Vg  =  vCog),  but  this  choice  will  lead  to  a  well-known 
problem,  termed  “sensitivity  to  mismatch"  by  Cox  [23].  Examples  of  this 

effect  will  be  given  in  the  next  section. 

For  the  present  study,  it  is  assumed  that  it  is  adequate  to  establish  a 
set  of  adaptive  weights  for  which  the  GIRCvg,  R,  v(ag))  exceeds  some 
minimal  operation  level,  such  as  (say)  10  dB.  Because  of  the  relatively  rapid 
convergence  of  this  performance  with  the  number  of  array  snapshots,  a 
reasonable  bound  on  the  adaptive  receiver  operating  characteristic,  may  be 
obtained  by  evaluating  the  GIR  for  Infinite  data.  In  the  examples  in  this 
report,  a  single  interferer  is  placed  A0  beamwidths  from  a  desired  signal  and 
the  locus  of  GIR  =  10  dB  for  AG  versus  desired  array  signal-to-nolse  ratio  is 
used  to  predict  limiting  performance. 
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III.  DESCRIPTION  OF  ALGORITHMS 


A.  Direction-Finding 

Many  practical  algorithms  that  simultaneously  estimate  the  directions  of 
several  emitters  generate  a  non-negative  function  called  a  DF  spectrum  from 
the  available  -  array  data.  The  domain  of  this  function  is  the  set  of  all 
possible  directions,  and  the  locations  of  its  maxima  (peaks)  correspond  to  the 
estimated  directions  of  arrival.  Most  algorithms  require  a  covariance 
estimate,  which  is  obtained  from  the  array  data.  The  covariance  estimate  is 

then  transformed  to  produce  the  spectral  estimate. 

For  the  purposes  of  this  discussion,  the  covariance  estimate  may  be  taken 
to  be  the  usual  sample  covariance  matrix  generated  from  K  snapshots  of  data 
{xk  I  k  =  1,  ...,  k},  i.e.. 


1  ^ 
I 


H 


R  ^  ^ 

k=l 


This  estimate  is  generally  positive  definite  provided  the  number  of  snapshots 
is  not  less  than  the  number  of  array  elements. 

Table  3.1  summarizes  the  type  of  transformation  used  by  several  of  the 
currently  most  popular  algorithms.  In  each  case,  the  algorithm  operates  on 
the  covariance  estimate  with  the  direction  vector  v(a)  to  generate  the 
spectral  estimate  for  the  direction  a.  As  can  be  seen,  all  of  the  DF  spectra 
in  table  3.1  require  the  computation  of  at  least  one  quadratic  form. 

For  the  case  of  a  uniform  linear  array,  any  quadratic  form 

g(  a)  =  v^  (  a)  G  v(  a) 

may  be  interpreted  as  a  "trig”  polynomial, .  provided  G  is  a  non-negative 
definite  Hermitian  matrix.  The  spectral  factorization  theorem  allows  us  to 
represent  g(a)  by  two  ordinary  polynomials,  one  with  its  roots  inside  the  unit 
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TABLE  3.1 


DIRECTION-FINDING  SPECTRA 


Adapted  Angular  Response 


v^(a)  R  ^  v(  g) 
v®(a)  r”^  v(a) 


Beam-Scan  Algorithm  (BSA) 

(  a)  R  v(  a) 

Maximum  Entropy  Method  (MEM) 


constant 
jv^(a)  R  ^ 

Maxinmm  Likllliood  Method  (MLM) 

_ 1 _ 

v®(a)  r"^  v(  a) 


u.is  the  first  column  of 
the  identity  matrix 


Multiple  Signal  Classification  (MUSIC) 


_ 1 _ 

v®(  a)  v(  a) 


certain  eigenvectors  of  R 
are  selected  for  the 
columns  of  E{j 


Thermal  Noise  Algorithm  (TNA) 

_ 1 _ 

v®(a)  R"^  v(a) 
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circle  and  the  other  with  corresponding  roots  outside  the  unit  circle. 
Moreover,,  if  z  is  a  root  of  one  of  these  polynomials,  then  1/z*  is  a  root  of 
the  other.  Either  of  these  two  polynomials  completely  characterizes  an 
"all-zero"  or  finite  impulse  response  (FIR)  filter  that  produces  a  random 
process  with  power  spectrum  g(a)  when  excited  with  white  noise.  A  process 
generated  in  this  manner  is  sometimes  called  a  moving  average  prbcess. 

The  beam  scan  algorithm  (BSA)  produces  a  spectral  estimate  consistent 
with  a  moving  average  process.  All  of  the  other  DF  spectra  in  Table  III  .A  are 
characterized  by  a  denominator  polynomial  except  AAR,  which  has  both  a 
denominator  and  a  numerator  polynomial.  The  AAR  spectrum  is  consistent  with 
an  autoregressive-moving  average  (ARMA)  process.  The  linear  filter  required 
to  synthesize  an  ARMA  process  has  both  poles  and  zeroes.  Each  of  the 
remaining  algorithms  (MEM,  MLM,  TNA,  and  MUSIC)  leads  to  a  spectral  estimate 
consistent  with  an  autoregressive  process.  As  one  might  expect,  an 
autoregressive  process  is  characterized  by  an  all-pole  filter. 

The  beam  scan  algorithm  actually  provides  the  best  possible  estimate  of 
the  direction  of  arrival  of  a  single  emitter  received  in  the  presence  of 
(spatially)  white  noise.  In  practice,  one  can  expect  the  beam  scan  method  to 
perform  adequately  so  long  as  the  emitters  are  all  well  Isolated.  Unfortu¬ 
nately,  this  approach  breaks  down  completely  when  two  (or  more)  emitters  are 

separated  by  less  than  an  array  beamwldth. 

The  beam  scan  algorithm  is  the  DF  equivalent  of  a  standard  technique  used 
in  classical  time— series  analysis.  In  many  applications,  a  single  record  of 
data  is  Fourier  transformed  to  obtain  a  rough  spectral  estimate  called  a  peri- 
odogram.  Averaging  many  periodograms  yields  a  smoothed  spectral  estimate. 
The  spectral  resolution  provided  by  this  approach  is  limited  by  the  length  of 
the  Individual  data  records.  Quite  analogously,  the  resolution  of  the  beam 
scan  algorithm  is  determined  by  the  length  of  the  array. 

The  limitations  of  the  traditional  approach  eventually  led  to  fundamental 
investigations  seeking  spectral  estimates  with  better  resolution.  Some  of  the 
more  Important  results  of  these  investigations  are  summarized  as  follows. 
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1.  Maximum  Entropy  Method 

In  1967,  John  Burg  [3]  shook  the  foundations  of  traditional  time-series 
analysis  with  his  assertion  that  conventional  spectral  estimation  techniques 
were  fundamentally  unsound.  Burg  was  upset  by  the  fact  that,  at  the  time,  all 
recognized  methods  for  computing  a  power  spectral  density  implicitly  truncated 
the  correlation  lags.  As  an  alternative.  Burg  proposed  his  now  famous 
••maximum  entropy"  spectral  estimate.  Although  it  may  not  have  been  widely 
recognized  at  the  time.  Burg  was  actually  advocating  that  spectral  estimates 
be  derived  within  the  framework  of  an  autoregressive  (AR)  model  for  the  time 

series  [4]. 

In  the  usual  time-series  setting,  maximum  entropy  spectral  estimates  are 
derived  from  a  linear  prediction  error  filter.  The  leading  coefficient  of 
this  filter  is  unity,  and  the  remaining  coefficients  are  chosen  to  minimize 
its  expected  output  power,  usually  referred  to  as  the  prediction  error.  The 
theoretical  basis  for  this  procedure  is  discussed  in  detail  in  Appendix  B  and 
the  references  therein.  The  desired  filter  coefficients  are  obtained  by  solv¬ 
ing  the  mixed  system  of  linear  equations 

5w  =  le  0  ...  Ol"^ 

for  the  prediction  error  e  and  the  unknown  elements  of 
w=[l?».*?]  » 

where  R  is  the  theoretical  covariance  (matrix)  for  an  arbitrary  snapshot.  The 
spectral  estimate  obtained  from  the  solution  to  (3.1)  can  be  written  as 

1 

S  -  —  ^ 

jv®(a)  R  Uj^j 
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where  ui  is  the  first  column  of  the  identity  matrix.  Strictly  speaking, 
this  expression  can  only  be  interpreted  as  a  maximum  entropy  spectral  density 
in  the  ideal  case  of  uncorrelated  emitters  and  a  linear  array  with  uniformly 
spaced,  omni-directional  elements.  Much  of  the  confusion  in  the  literature 
concerning  the  maximum  entropy  approach  (see,  for  example,  [5])  can  be 
attributed  to  the  prevailing  uncertainty  regarding  the  proper  choice  of  a 
covariance  estimate.  The  elegance  and  efficiency  of  the  well-known  Levinson 
recursion  [6]  has  prompted  many  researchers  to  force  a  Toeplitz  structure  on 
the  covariance  estimate.  The  standard  Yule-Walker  method  [7]  leads  to  an 
estimate  that  is  positive  definite  but  badly  biased,  at  least  for  DF 
applications.  The  bias  is  easily  removed,  but  only  at  the  expense  of 

destroying  (with  some  non-zero  probability)  the  desired  positive  definite 
property.  Neither  of  these  approaches  is  recommended  when  high  resolution  is 

important • 

Fortunately,  the  intrinsic  (temporal)  smoothing  provided  by  a  sample  co- 
variance  matrix  is  quite  adequate  for  most  direction-finding  applications.  Of 
course,  sample  covariance  matrices  are  never  Toeplitz  (except  by  accident) 
and,  if  employed,  one  must  then  solve  the  linear  prediction  (3.1)  without  the 
help  of  a  truly  fast  algorithm.  However,  in  most  DF  applications,  the  data 
records  (i.e.,  snapshots)  are  short  and  computational  efficiency  is  relatively 
important.  In  the  (rare)  situations  where  the  sample  covariance  approach  is 
unsatisfactory,  the  problems  encountered  with  the  Yule-Walker  approach  can  be 
circumvented  by  Burg's  ingenious  scheme  for  estimating  reflection  coefficients 
directly  from  the  data. 

The  standard  time-series  implementation  of  Burg's  technique  [8]  is  quite 
efficient,  and  a  Burg  filter  always  has  the  minimum  phase  property.  This 
property  guarantees  a  stable  inverse  but  is  seldom  crucial  except  perhaps  when 
one  wishes  to  synthesize  the  input  process.  Moreover,  the  extended  Burg  tech¬ 
nique-  (9]  for  processing  multiple  snapshots  is  not  significantly  faster  than 
standard  (e.g.,  Cholesky  type)  algorithms  for  solving  the  prediction  (3.1)  via 
the  sample  covariance  method  [10]. 
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2*  Maximum  Likelihood  Method 

In  [11]  Lacoss  discusses  the  maximum  entropy  method  (liEM)  and  another 
high-resolution  spectral  estimation  algorithm  attributed  to  J.  Capon  [12] 
called  the  maximum  likelihood  method  (MLM) .  The  basic  idea  behind  the  latter 
approach  is  simple.  Weights  for  the  array  elements  are  chosen  which  insure 
unit  gain  in  a  given  direction  o  while  simultaneously  minimizing  the  array 
output  power.  Under  these  conditions,  the  output  of  the  adapted  array  pro¬ 
vides  an  unbiased,  minimum  variance  estimate  of  the  (desired)  signal  arriving 
from  the  specified  direction.  When  the  interference  is  Gaussian,  the  power 
out  of  the  adapted  array  is  a  maximum  likelihood,  estimate  of  the  power  re- 
celved  from  the  direction  ou 

In  mathematical  terms,  the  power  out  of  the  array  can  be  expressed  as 


w  represents  the  array  weights.  Constraining  the  gain  of  the  array  to 
be  unity  in  the  direction  a  is  achieved  by  demanding  that  w  satisfy 

w^ v(a)  ®  1  .  (3.3) 

Minimizing  (3.2)  subject  to  (3.3)  is  easily  accomplished  using  the  method  of 
LaGrange  multipliers.  The  optimum  weights  for  this  problem  maximize  the 
output  signal-to-interference  ratio  (SIR)  and  are  found  by  solving 


Rw(  a)  =  X  v(  a)  (3.4) 

where  the  Lagrangian  X  is  chosen  to  satisfy  (3.3).  Pre— multiplying  this 
equation  by  the  Hermitian  transpose  of  w(a),  we  find  that  the  MLM  power 
estimate  and  the  Lagrangian  are  numerically  the  same.  Replacing  X  with  P  in 
(3.4),  solving  for  w(a),  and  substituting  the  result  in  (3.2)  leads  to  the  MLM 
power  estimate 


20 


1 


MLM 


v^(a)  R  ^v(a) 


3.  Adapted  Angular  Response 

The  adapted  array  response  (AAR)  algorithm  suggested  by  Borgiottia  and 
Kaplan  [13]  can  be  interpreted  as  a  variation  on  the  MLM  theme.  As  mentioned 
above,  the  MLM  array  weights  maximize  the  output  SIR.  This  fact  remains  true 
for  any  choice  of  the  Lagrange  multiplier.  The  AAR  spectral  estimate  is  gen¬ 
erated  by  scaling  the  MLM  weights  so  that  the  sum  of  their  squared  magnitudes 
is  some  fixed  value,  e.g.. 


w(a) 


=  1 


(3.5) 


This  modification  leads  to  a  DF  method  with  the  desirable  property  that  the 
effect  of  white  noise  on  the  spectral  estimate  is  (on  the  average)  the  same  in 
every  direction.  Choosing  the  Lagrangian  X  to  satisfy  (3.5)  instead  of  (3.3) 
leads  to  the  AAR  power  estimate 


v^(a)  R  ^v(a) 


v^a)R-2v(a) 


4.  Thermal  Noise  Algorithm 

As  was  mentioned  above,  the  AAR  power  estimate  has  the  same  mathematical 
structure  as  an  ARMA  power  spectral  density.  Generally  speaking,  an  ARMA 
process  can  be  represented  by  the  cascade  of  an  all-zero  (FIR)  filter  and  an 
all-pole  filter.  In  [14],  the  point  was  made  that  a  signal  consisting  of  sin¬ 
usoids  (i.e.,  plane  waves)  and  additive  white  noise  satisfies  an  ARMA— li*e 


21 


difference  equation.  However,  the  sinusoids  in  noise  process  is  an  extremely 
pathological  case  where  the  poles  and  zeroes  lie  on  the  unit  circle  and  cancel 
exactly.  Thus,  the  frequency  response  of  the  cascade  is  perfectly  flat!  The 
sinusoids  arise  from  the  transient  response  of  the  critically  stable  all-pole 
filter.  This  strongly  suggests  that  the  denominator  (l.e.,  the  AR  part)  of  an 
ARMA  spectral  estimate  suffices  to  determine  the  frequencies  of  the  sinu¬ 
soids.  Although  the  underlying  reasons  that  motivated  Gabriel  [15]  were  un¬ 
doubtedly  somewhat  different,  the  thermal  noise  algorithm  (TNA)  uses  only  the 
denominator  of  the  AAR  spectrum,  l.e., 

‘  .“(a)  v(a) 


5.  Multiple  signal  Classification 

The  Multiple  Signal  Classification  (MUSIC)  approach  to  direction-finding 
was  first  described  in  [16].  The  theoretical  framework  behind  the  MUSIC  algo¬ 
rithm  [17]  Is  quite  general  and  substantially  extends  the  pioneering  harmonic 
retrieval  method  of  Pisarenko  1 18 ] . 

The  underlying  assumption  behind  the  MUSIC  algorithm  is  that  the  numoer 
of  emitters  seen  by  the  receiver  is  less  than  the  number  of  antenna  elements. 
Under  this  condition,  the  covariance  matrix  of  the  received  signal 


S  =  e{ss®} 

is  singular.  Referring  to  (3.1),  we  observe  that  the  null  vectors  of  S 
theoretically  provide  a  perfect  mechanism  for  spatially  extrapolating  (i.e., 
"predicting")  the  received  signal.  Thus,  one  intuitively  expects  that  the 
emitter  directions  could  somehow  be  extracted  from  the  null  space  of  S. 
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The  space  spanned  by  the  columns  of  S  is  referred  to  as  the  signal 
space.  Except  in  certain  pathological  cases  (i.e.,  perfectly  correlated  sig¬ 
nals),  the  direction  vectors  of  the  emitters  will  always  He  in  the  signal 
space.  On  the  other  hand,  an  arbitrary  direction  vector  will  generally  have  a 
component  in  the  (complementary)  null  space.  Thus,  a  simple  test  based  on  the 
distance  to  the  signal  space  determines  whether  or  not  an  arbitrarily  chosen 
direction  Is  a  possible  emitter  direction. 

Given  an  orthonormal  basis  for  the  null  space,  the  Euclidian  distance 
from  an  arbitrary  vector  x  to  the  signal  space  can  be  easily  calculated.  The 
projection  of  x  into  the  null  space  can  be  written  as 

where  the  columns  of  E^  are  the  (orthonormal)  basis  vectors  for  the  null 
space.  Consequently,  the  distance  from  x  to  the  signal  space  is 

Kl  ■  ’ 

since 


by  construction* 

The  MUSIC  (pseudo)  spectrum  is  defined  to  be  the  inverse  of  the  squared 
distance  from  an  arbitrary  direction  vector  to  the  signal  space,  i.e.. 


l^MUSIC 


_ 1 _ 

v®(  o)  Ej^®  v(  a) 
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In  theory,  the  MUSIC  spectrum  becomes  infinite  at  the  true  directions.  How¬ 
ever,  small  errors  in  the  null  vectors  or  the  direction  vectors  will  almost 
surely  prevent  the  denominator  from  vanishing  entirely.  Therefore,  in  prac¬ 
tice,  the  direction  vectors  that  lie  closer  to  the  signal  space  than  their  Im¬ 
mediate  neighbors  determine  the  MUSIC  direction  estimates. 

In  an  operational  system,  the  MUSIC  spectrum  must  be  derived  from  a 
sample  covariance  matrix.  The  first  step  In  the  MUSIC  algorithm  Is  to  deter¬ 
mine  the  standard  eigenvalue  (spectral)  decomposition  [19]  of  the  sample  co- 
variance  matrix,  l»e., 

R  =  EDE 

where  E  Is  a  unitary  matrix  and  D  Is  a  (positive)  diagonal  matrix.  Post- 
multiplying  by  E,  we  Immediately  observe  that  the  columns  of  E  are  eigenvec¬ 
tors  of  the  sample  covariance  matrix,  and  the  (diagonal)  elements  of  D  are  the 
corresponding  eigenvalues.  Without  loss  of  generality,  we  may  assume  that  the 
eigenvalues  have  been  arranged  in  descending  numerical  order. 

The  signal  space  is  specified  by  a  simple  partition  of  E,  i.e., 

E  =  [Eg  I  E^l 

At  this  point,  the  critical  issue  is  the  proper  dimension  of  the  signal 
space.  If  the  number  of  emitters  present  is  known  apriori  (never  the  case  in 
practice),  then  the  correct  partition  to  employ  Is  obvious.  Given  J  emitters. 
Eg  consists  of  the  first  J  columns  of  E.  However,  when  the  number  of  emitt¬ 
ers  is  unknown,  a  choice  (for  J)  must  be  made  based  on  the  available  data. 
Several  methods  for  determining  J  from  the  eigenvalues  have  been  examined. 
These  algorithms  perform  a  sequence  of  likelihood  ratio  (hypothesis)  tests  and 
select  the  smallest  value  of  J  that  is  statistically  consistent  with  the  em¬ 
pirical  eigenvalue  distribution.  Simulations  results  and  references  for  the 
various  likelihood  ratio  tests  studied  thus  far  are  presented  in  Appendix  C. 
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6.  Rooting  Methods  for  Linear  Arrays 
The  root-finding  method  discussed  here  applies  to  all-pole  (autoregres¬ 
sive)  spectra  of  the  general  form 


S(a)  = 


1 

v^(a)  P  v(a) 


where  P  is  a  non-negative  definite  KyU  Hermitian  matrix  that  generally  depends 
on  the  sample  covariance  matrix  in  a  nonlinear  manner.  For  example,  the  lOK 
spectrum  is  obtained  by  choosing  P  to  be  the  Inverse  of  the  (sample) 
covariance  matrix.  For  the  purposes  of  this  report,  we  may  restrict  our  at¬ 
tention  to  uniform  linear  arrays.  In  this  special  case,  the  direction  vectors 
are  specified  in  Section  IV. A.  However,  the  basic  approach  described  here  is 
easily  extended  to  thinned  (uniform)  arrays  with  missing  elements. 

A  direct  calculation  shows  that  the  inverse  of  the  (autoregressive) 

spectrum  defined  above  can  be  written  as 
M-1 

^  p  exp  {-i  2irai  » 

m  =  “M+1 

where  the  mth  coefficient  is  obtained  by  summing  the  elements  on  the  mth  diag- 
onal  of  P,  i.e*. 


P 


m 


M 


I 

k-l=m 


^kl 


(3.6) 


By  introducing  the  change  of  variables 


z 


exp  f  i  2it  ^a} 


(3.7) 
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the  spectrum  can  be  expressed  in  terms  of  the  polynomial 


.  M— 1  v  “in 

P(z)  =  z  I  p  z 

m  =  -M+1 

Specifically,  we  have 

S  =  jpfexp  {-i2TrCa}}j  ^ 

Clearly,  the  polynomial  P(2)  contains  all  the  information  In  the  spectrum.  In 
fact,  spectral  peaks  are  a  natural  consequence  of  polynomial  roots  that  lie 
near  the  unit  circle.  Hence,  we  are  led  to  explore  rootfinding  as  an  altema“ 
tive  to  searching  for  spectral  peaks. 

The  Hermltian  property  of  P  insures  that  the  polynomial  coefficients  in 

(3.6)  satisfy  the  (harmonic)  relationship 

* 

p  =  p 

-m  m 

Using  this  fact,  it  can  be  verified  that: 

★ 

If  z  is  a  root  of  P(2),  so  is  1/z  . 

m  ra 

In  practice,  exactly  half  the  roots  of  P(z)  will  lie  inside  the  unit  circle. 
Direction  (cosine)  estimates  can  be  derived  from  the  roots  by  referring  to 

(3.7) ,  i.e., 

-  _  ^m 

**m  2ir  5 
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When  the  separation  between  adjacent  elements  is  less  than  1/2  wavelength, 
(i.e.-,  ?  <  1/2),  the  magnitude  of  an  estimated  direction  cosine  may  be  larger 
than  unity.  In  this  case,  the  estimate  lies  outside  "visible"  space  and 
should  be  discarded.  Naturally,  other  criteria  can  also  be  used  to  eliminate 
spurious  estimates.  For  example,  roots  that  lie  close  to  the  origin  are  pre¬ 
sumably  of  little  interest,  since  they  do  not  correspond  to  (significant) 
spectral  peaks.  In  the  case  of  the  root  version  of  the  MUSIC  algorithm,  only 
the  J  roots  closest  to  the  unit  circle  are  considered,  where  J  is  the  (esti¬ 
mated)  dimension  of  the  signal  space. 

B.  Power  Estimation 

It  is  well-known  that  sampling  of  a  time  function  introduces  ambiguities 
in  the  resulting  spectrum.  Similarly,  sampling  the  received  wavefront  on  an 
aperture  produces  ambiguities  in  the  resulting  estimates  of  angle  of 
arrival.  The  most  familiar  example  is  that  of  grating  lobes  which  are 
produced  when  the  elements  of  an  array  are  spaced  farther  than  X/2  apart. 
The  ambiguity  introduced  by  a  grating  lobe  is  fundamental;  there  is  no  way  to 
determine  from  which  of  the  two  (or  more)  directions  signals  are  arriving. 

The  existence  of  ambiguities  is  a  consequence  of  the  linear  dependence 
of  direction  vectors.  For  example,  suppose  there  exist  direction  vectors 
,  y2  >  ^  such  that 

av  ^  +  by^2  ® 

for  some  set  of  complex  numbers  (a,b,c),  and  that  the  signal  is  composed  of 
sources  from  directions  complex  Gaussian  amplitudes  A,  B»  If 

the  correct  directions  are  selected,  the  resulting  signal-in-space  covariance 
matrix  is 


r*  A  “ 

0 

[b_ 

(A  B  )  = 

0 

1 _ 
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The  signal  can  be  written  equally  well  as 


jT  =  (B  -  aA)v2  -  bAv^2 

However,  if  directions  V3  are  selected,  the  resulting  signal-in-space 

covariance  is 


which  is  not  diagonal.  The  same  is  true  if  directions  vj,  3,  or  vj, 
V9 ,  vq  are  chosen. 

Adaptive  algorithms  for  estimating  the  direction  of  arrival  of 
incoherent  signals  all  operate  on  the  principle  of  adaptive  cancellation. 
They  select  a  set  (or  sets)  of  weights  to  apply  to  the  array  data  so  as  to 
minimize  the  output  power  from  the  array.  Some  sort  of  constraint  must  be 
Imposed  on  the  weights  to  keep  them  from  going  to  zero;  the  algorithms  differ 
mainly  in  their  choice  of  this  constraint. 

It  is  intuitively  clear  that  the  adapted  array  pattern  (i.e.,  the  array 
pattern  using  the  adaptive  weights)  must  have  minimal  gain  in  each  of  the 
source  directions.  These  minima  serve  as  estimates  of  the  source 
directions.  If  32,  33  are  linearly  dependent  directions,  and 

sources  are  present  in  directions  vj,  32.  the  adapted  pattern  will  have 
minima  in  all  three  directions.  To  determine  the  true  source  distribution,  a 
power  estimation  algorithm  which  can  handle  linearly  dependent  direction 
vectors  is  needed. 

The  usual  method  of  power  estimation  requires  that  the  candidate 
direction  vectors  be  linearly  independent.  Then,  knowing  the  noise  power 
0^,  one  can  operate  on  the  tirue  covariance  matrix 
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R  =  +  VPV® 


and  obtain 

P  =  (v®V)"^  (R  -  a^I)  V  (V®V)"^  (3.8) 

Of  course,  R  is  not  known,  so  that  sample  covariance  matrix 


H 

r  r 
— n-n 


(3.9) 


is  used  instead.  If  the  direction  vectors  are  linearly  dependent, 
(vHv)“l  does  not  exist,  and  the  method  fails. 

There  are  a  variety  of  ways  to  avoid  this  problem.  Several  are 
described  in  a  recent  paper  by  d’Assumpcao  [20].  The  technique  discussed  here 
is  called  least  squares  power  estimation. 

The  idea  of  least  squares  power  estimation  is  to  find  positive  numbers 
^1»  *  •  noise  power  is  not  known)  which  minimize 


E  =  ns  -  I  -  VPvHp,^ 


(3.10) 


=  Tr(S  -  -  VPV®)“  (S  -  -  VPV®) 

=  11.  |(S  -  A  .  ,pvH)^j2 


Assume  first  that  a  is  knowtiv 
set  of  equations^ 


Setting  =  0  for  £==1,L,  we  obtain  the 


1 

Note  that 


(VPV^) 


H 

V  V  • 
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(S  -  ch)v^  =  VPV^ 


Z=l,L 


L  .  E  2 

’  Vl 


These  equations  may  be  written  in  matrix  form  as 


OP  =  b  -  0  c 


(3.11) 


where 


=  Zi  S 


(3.12) 


‘'A  “ 


(3.13) 


«  H  2 


(3.14) 


The  matrix  Q  can  be  written  as  a  Hadamard  product.  The  Hadamard  product  of 
two  matrices  is  defined  to  be 


(A  □  B)^^  ®k£\£ 


(3.15) 


With  this  notation  we  can  write 


0  =  V®v  □  (v“v)’ 


(3.16) 
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(3.17) 


2 

where  *  denotes  the  complex  conjugate. 

The  solution  of  (3.11)  is  clearly 

^  =  o”^(h  -  a^£) 

In  many  cases,  0~1  exists  even  when  (V^V)”!  does  not. 

When  is  unknown,  we  first  solve 

il.  =  -2  Tr  (S  -  a^I  -  VPV^)  =  0 

with  the  result 

=  4  TR(S-VPV®) 

n 

where  M  is  the  number  of  array  elements  (the  dimension  of  S) 
into  (3.10)  gives 

E  =  Tr(S  -  0^1  -  VPV^)®  (S  -  vpv") 

=  Tr(S  -  VPV®)"  (S  -  vpv")  -  4  Tr^(S  -  VPV®) 


Setting  the  derivatives  equal  to  zero  yields 


u 

k+1 


^  (  V^  V 


TrS 

M 


^ince  V®V  is  Herraitian,  we  could  equally  well  write 
0  =  #V  X  (VHV)T. 


(3.18) 

.  Substitution 
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wliich  C3.n  be  writteri  in.  tbe  form 


f  r\  ^  T) 

(0  -  M  £  £  )  I 


(3.19) 


The  solution  is,  of  course 


P 


/n  ^  (\\ 

(0  -  ^  £  £  ) 


(3.20) 


If  0  has  an  inverse,  so  does  0  “  ^  £  £  • 

Unfortunately,  (3.17),  (3.18),  and  (3.20)  do  not  guarantee  positive 
estimates  for  the  signal  and  noise  powers.  This  reflects  the  fact  that  in  the 
minimization,  no  constraint  on  the  parameter  values  was  imposed.  If  the 
minimization  were  done  numerically  hy  a  nonlinear  program,  the  constraints 
could  easily  be  included.  However,  this  would  be  computationally  much  more 
expensive  than  the  simple  unconstrained  solution. 

The  least  squares  power  estimation  technique  eliminates  spurious 
spectral  peaks  by  producing  small  (perhaps  even  negative)  estimates  of  the 
power  associated  with  them;  they  can  then  be  eliminated  by  a  simple  threshold 
test.  The  results  of  a  Monte  Carlo  computer  simulation  comparing  the  direct 
method  ((3.8),  retaining  only  the  diagonal  terms)  and  the  least  squares  method 
(3.17)  are  presented  in  Appendix  D. 


C.  Weight  Design  Procedures  for  Adaptive  Listening 

The  problem  of  designing  a  set  of  optimal  weights  for  an  adaptive  lis¬ 
tening  array  has  been  considered.  This  problem  has  been  treated  extensively 
in  the  past  [21-28].  In  these  previous  analyses,  it  has  been  assumed  that  the 
covariance  matrix  of  the  interference,  or  noise,  can  be  obtained,  or  is 
possibly  known,  in  the  absence  of  the  desired  signal.  However,  in  the  present 
work,  it  is  assumed  that  the  covariance  matrix  of  the  signal  plus  noise  is  the 
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quantity  which  can  be  obtained.  This  assumption  leads  to  great  complications, 
as  has  been  observed  by  Cox  [23].  In  particular,  if  an  attempt 
is  made  to  use  this  covariance  matrix  to  design  a  set  of  weights,  using  the 
procedure  which  maximizes  the  signal-to-interference  ratio  (SIR),  then  the 
resultant  processor  will  be  extremely  sensitive  to  the  mismatch  problem, 
i.e.,  the  assumption  about  the  angle  of  arrival  of  the  desired  signal.  This 
particular  design  procedure  is  termed  the  Measured  Covariance  Method,  and  the 

weights  are  designed  as 


where  Rx  is  the  covariance  matrix  for  the  desired  signal  plus  interference, 
discussed  in  Section  II. B.,  and  V  is  a  steering  vector  which  points  in  the 

estimated  direction  of  the  desired  signal* 

A  number  of  other  design  procedures  were  considered.  The  first  method  is 

termed  the  Model  Covariance  Method  and  the  weights  are  computed  as 


where  J  is  the  number  of  emitters,  is  a  steering  vector  which  points  in 
the  estimated  direction  of  the  k-th  emitter,  <P-  is  the  incoherent  noise 
power,  I  is  the  L  X  L  identity  matrix.  It  is  assumed  that  the  desired  emitter 
corresponds  to  the  k-th  direction.  In  practice,  the  steering  vectors  would  be 
obtained  from  a  direction-finding  algorithm. 

An  example  of  the  problem  of  sensitivity  to  mismatch  is  shown  in  Figs. 
3.1a-b.  A  linear  array  of  10  elements  is  considered  with  an  inter-element 
spacing  of  0.5  wavelengths.  The  power  levels  of  the  desired  and  undesired 
emitters,  relative  to  the  background  incoherent  noise,  is  26  dB  and  40  dB, 
respectively,  per  array  element. 


33 


EMITTER  SEPARATION  (Beamwidth) 


Fig.  3.1.  Copy  sensitivity  to  signal  direction-of-arrival  estimation  errors. 

a)  Typical  beam  patterns  with  emitter  separation. 

b)  Typical  array  output  performance  vs.  emitter  separation. 
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In  Tig.  3.1a  the  beam  patterns  are  depicted  for  the  Measured  Covariance 
Method  and  the  Model  Covariance  method.  The  directions  of  arrival  of  the 
desired  emitter,  and  interference,  are  assumed  to  be  0  and  0.5  beamwidths, 
respectively.  The  pattern  for  the  sum  beam,  corresponding  to  a  simple  phase 
steering  and  then  summing  of  the  antenna  element  outputs,  serves  to  define  the 
natural  beamwidth  of  the  array  as  determined  from  the  dimension  of  the  array 
in  wavelengths,  which  was  discussed  in  Section  II. B.  The  pattern  for  the 
Measured  Covariance  Method,  when  there  is  no  error  in  the  estimated  direction 
of  arrival  of  the  desired  emitter,  produces  a  null  in  the  direction  of  the 
interference,  but  has  a  relative  array  gain  of  about  0  dB  for  the  desired 
signal.  However,  if  an  error  of  0.05  beamwidths  is  made,  the  pattern  for  the 
Measured  Covariance  Method  produces  a  null  in  the  direction  of  the  desired 
signal.  The  reason  for  this  is  that  the  design  procedure  treats  the  desired 
signal  as  interference,  since  its  direction  of  arrival  has  not  been  specified 
precisely.  The  pattern  for  the  Model  Covariance  Method  under  these  conditions 
does  not  exhibit  such  behavior  as  seen  in  Fig.  3.1a.  This  pattern  produces  a 
null  in  the  direction  of  the  interference,  but  passes  the  desired  signal.  In 
essence,  this  behavior  of  the  pattern  occurs  since  the  effect  of  the  desired 
signal  on  the  covariance  matrix  used  in  the  design  procedure  has  been 
eliminated.  Thus,  the  Model  Covariance  Method  is  relatively  insensitive  to 


signal  DOA  errors. 

The  effect  of  emitter  separation  upon  the  sensitivity  to  mismatch  problem 
is  shown  in  Fig.  3.1b  which  depicts  the  output  signal-to-interference  ratio 
(SIR)  in  dB,  vs.  the  emitter  separation  in  beamwidths.  The  method  for 
computing  the  SIR  was  discussed  in  Section  II. B.  The  result  for  the  sum  beam 
shown  in  this  figure  indicates  that  the  performance  obtained  with  this  method 
is  not  impressive  when  the  magnitude  of  the  emitter  separation  is  less  than  1 
beamwidth.  However,  optimal  performance  is  obtained  for  the  Measured 
Covariance  Method  when  there  is  no  error  made  in  estimating  the  direction  of 
arrival  of  the  desired  signal.  If  the  emitter  separation  is  0.5  beamwidths, 
for  example,  then  it  is  possible  to  obtain  a  depth  of  null  of  about  45  dB. 


That  is,  the  desired  signal  level  can  be  raised  from  14  dB  below  the 
Interference  to  about  31  dB  above  the  interference.  Unfortunately,  if  an 
error  of  0.05  beamwidths  is  incurred,  then  the  performance  of  this  method 
deteriorates  seriously,  as  shown  in  Fig.  3.1b.  However,  the  impressive 
performance  is  maintained  by  the  Model  Covariance  Method  under  these 
conditions,  as  depicted  in  Fig.  lb.  These  results  indicate,  once  again,  that 
the  Model  Covariance  Method  is  not  sensitive  to  the  error  made  in  estimating 

the  direction  of  arrival  of  the  desired  signal. 

Another  weight  design  procedure  which  has  been  considered  is  the  Projec¬ 
tion  Nulling  Method.  In  this  method  the  complex  weights  are  computed  as 

Wp  =  (I  -  V(V^V)“^V®)  V^. 


where 


V  =  [Vj 


It  has  been  shown  ([8],  p.  141)  that  if  the  slgnal-to-noise  ratio  is  large, 
then  the  Projection  Nulling  Method  is  equivalent  to  the  Model  Covariance 
Method.  This  theoretical  result  has  been  confirmed  by  computer  simulations 
in  the  present  study.  As  a  consequence,  the  Projection  Nulling  Method  has  not 
been  discussed  in  detail.  As  a  final  comment,  it  should  be  noted  that  since 
the  Projection  Nulling  Method  is  equivalent  to  the  Model  Covariance  Method, 
the  former  method  is  also  not  sensitive  to  the  error  incurred  in  measuring  the 
direction  of  arrival  of  the  desired  signal. 
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IV.  DESCRIPTION  OF  THE  MONTE  CARLO  EXPERIMENTS 


A.  Direction-Finding  Comparison 

A  primary  motivating  factor  behind  this  report  is  the  desire  to 
communicate  the  results  of  an  extensive  simulation  of  a  number  of  potentially 
interesting  direction-finding  algorithms.  These  quantitative  performance 
comparisons  are  presented,  in  their  entirety,  in  Appendix.  F.  A  conceptual 
flow  diagram  of  the  simulation  experiments  is  shown  in  Fig.  4.1. 

The  basic  scenario  for  all  of  the  experiments  consisted  of  two 
independent  emitters  and  a  uniform  linear  array  with  ten  omni-directional 
elements.  Adjacent  elements  were  separated  by  1/2  wavelength.  The  only 
source  of  error  was  additive  thermal  noise  at  each  of  the  array  elements. 
The  signals  and  noise  were  modelled  as  complex  Gaussian  variables  as 
discussed  in  Section  II.A.  Sample  covariance  matrices  were  generated  using 
the  Wishart  technique  described  in  Appendix  E.  Since  our  antenna  model 
exhibits  perfect  symmetry  with  respect  to  its  geometric  (phase)  center,  the 
sample  covariance  matrix  was  first  processed  using  the  forward/backward 
averaging  technique  described  below.  Based  on  the  resulting  covariance 
estimate,  tentative  direction  of  arrival  estimates  were  generated  using  one 
of  the  algorithms  discussed  in  Section  III .A.  The  array  slgnal-to-noise 
ratio  (SNR)  associated  with  each  candidate  direction  was  estimated  using  the 
direct  method  discussed  in  Section  III.B.  All  candidate  directions  with  an 
SNR  estimate  less  than  0  dB  were  discarded.  The  remaining  direction 

estimates  were  compared  with  the  true  directions.  Those  direction  estimates 
resulting  in  the  smallest  total  (absolute)  direction  error  were  assigned  to 
the  true  directions  subject  to  the  following  provisions: 

1 .  A  single  estimate  could  not  be  assigned  to  both 
emitters. 

2.  Assignments  with  an  absolute  direction  error  larger  than 
a  beamwidth  were  not  permitted. 

3.  An  assignment  was  disallowed  if  the  power  estimate  was 
10  dB  larger  than  the  true  power. 
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Fig,  4.1.  Process  for  simulating  superresolution. 
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All  unassigned  estimates  were  counted  as  false  alarms.  An  unassigned  emitter 
was  declared  to  be  a  miss.  The  emitters  were  said  to  be  resolved  if  neither 
was  missed.  The  probability  of  resolution  and  the  false  alarm  rate  were 
calculated  based  on  100  Monte  Carlo  trials.  Second  order  statistics  on  the 
DF  errors  were  also  accumulated,  conditioned  upon  successfully  resolving  the 

two  emitters. 

During  the  course  of  the  simulation  study,  the  fundamental  parameters  of 
the  system  were  varied  in  order  to  provide  a  better  perspective  on  the 
relative  performance  of  the  algorithms.  The  sensitivity  of  the  algorithms  to 
thermal  noise  was  tested  by  changing  the  number  of  snapshots  used  to 
construct  the  sample  covariance.  Three  cases  were  considered,  namely  10, 
100,  and  1000  snapshots  (looks).  To  assess  the  resolution  limits  of  the  DF 
algorithms,  emitter  separations  of  0.1,  0.2,  and  0.4  beamwidths  were 
considered.  Initially,  the  relative  power  of  the  two  emitters  was  0  dB. 
Subsequently,  an  identical  series  of  experiments  was  conducted  with  the 
relative  power  set  at  -10  dB.  In  both  cases,  the  array  SNR  of  the  desired 
(weaker)  signal  was  varied  from  10  to  50  dB  in  5  dB  steps.  These  parameter 
variations  are  summarized  in  Table  4.1. 


Table  4.1 


PARAMETER 

VARIATIONS 

SIR  : 

0,  -10  dB 

Emitter  Separation  : 

0.1,  0.2,  0.4  beamwidths 

Number  of  Looks  : 

10,  100,  1000 

Array  SNR  : 

10,  15,  ...,  45,  50  dB 
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The  relative  power  of  the  desired  signal  to  the  interfering  signal  is  called 
the  signal-to-interf erence  ratio  (SIR).  The  emitter  separation  is  specified 
in  terms  of  the  available  beamwidth  of  the  array,  defined  to  be  the  inverse 
of  the  length  of  the  array  expressed  in  wavelengths.  The  desired  signal  was 
always  located  broadside  to  the  array  (i.e.,  o  =  0),  and  the  beamwidth  was 
0,2  radians  or  11.5  degrees.  The  number  of  independent  (uncorrelated) 
snapshots  used  to  construct  the  sample  covariance  matrix,  denoted  by  K,  is 
referred  to  as  the  number  of  "looks".  If  the  receiver  IF  bandwidth  is  B  and 
the  data  collection  interval  is  T,  then  the  number  of  looks  is  limited  by  the 
(time“bandwldth)  product  BT,  i.e., 

K  £  BT  . 

The  array  SNR  is  the  total  energy  (in  joules)  received  from  an  emitter  in  a 
single  snapshot  divided  by  the  thermal  noise  level  (in  watts/Hz  =  joules). 
Note  that  the  array  SNR  includes  the  signal  gain  available  from  the  antenna. 
For  the  ideal  ten  element  array  considered  here,  the  array  SNR  is  10  dB 

greater  than  the  SNR  on  a  single  antenna  element. 

Calculations  of  the  array  SNR  are  facilitated  by  scaling  the  direction 
vector  v(a)  to  have  unit  norm.  Moreover,  it  is  usually  convenient  to  choose 
the  phase  reference  (center)  to  be  at  the  geometric  center  of  the  array. 
Thus,  the  mth  element  of  the  normalized  direction  vector  for  an  ideal  array 
is 


V  (a)  = 


M 


1/2 


exp  {i  2'ir(m  —  M/2)  5®}  >  m  -  0,  1, 


M-1 


(4.1) 


where  M  is  the  number  of  array  elements  and  C  (-  d/X)  is  Che  element 
separation  in  wavelengths.  One  may  easily  verify  that  Che  direction  vector 
constructed  from  (4.1)  has  the  properties 


2 


1 
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and 


v(a)  =  Jv  (a)  ,  (^*2) 

where  J  is  the  (usual)  exchange  matrix  that  reverses  the  order  of  the 

elements  of  a  vector  (see  Appendix  B). 

A  vector  that  satisfies  the  relationship  in  (4.2)  is  said  to  be 
harmonic.  A  simple  calculation  reveals  that  every  symmetric  array  with 
isotropic  (i.e.,  identical)  elements  has  harmonic  direction  vectors .  In 
turn,  the  covariance  matrix  of  the  received  signal,  s  =  Vp  (see  Section 
II. A),  satisfies  a  similar  harmonic  property,  e.g., 

S  =  e{ss^} 

=  VPV® 


*  T 

=  JV  PV  J 

* 

=  JS  J 

where  P  is  a  real  diagonal  matrix  consisting  of  the  emitter  powers  that  would 
be  measured  at  the  array  phase  center.  When  a  Gaussian  signal  with  a 
harmonic  (persymmetrlc)  covariance  matrix  is  received  in  Gaussian  white 
noise,  it  can  be  shown  [29]  that  the  forward/backward  averaged  sample 
covariance  matrix,  i-e., 

.  (*•3) 
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Fig.  4.2.  Typical  output  from  a  superresolution  experiment 


is  a  sufficient  statistic  for  estimating  any  parameter  of  the  (true) 
covariance  matrix.  Note  that  this  (harmonic)  covariance  estimate  is  obtained 
by  taking  the  arithmetic  average  of  the  usual  sample  covariance  matrix  and 
the  (sample)  covariance  matrix  constructed  from  reversed  and  conjugated 
data.  All  of  the  algorithms  compared  in  this  report  derived  their  DF 

estimates  from  the  covariance  estimate  in  (4.3). 

Examples  of  the  statistical  data  produced  by  the  simulation  are  shown  in 
Fig.  4.2.  Each  data  point  was  calculated  on  the  basis  of  100  Monte  Carlo 
trials.  Smooth  curves  were  fitted  to  the  simulation  data  with  a  standard 
algorithm  provided  in  the  software  package  used  to  plot  the  results.  For 
purely  aesthetic  reasons,  the  actual  data  points  are  normally  suppressed. 

Perhaps  the  most  important  DF  statistic  is  the  probability  of  resolving 
two  emitters.  Generally  speaking,  most  DF  algorithms  cannot  reliably  resolve 
closely  spaced  emitters  at  very  low  signal-to-nolse  ratios.  However,  as 
indicated  by  the  example  in  Fig.  4.2,  an  algorithm’s  ability  to  resolve 
emitters  improves  rapidly  once  the  SNR  exceeds  some  critical  threshold 
level.  Below  threshold,  most  algorithms  erroneously  report  the  presence  of  a 
single  emitter  at  the  "centroid"  of  the  two  emitters.  At  the  two  extremes, 
relatively  few  spurious  estimates  are  generated.  However,  the  transition 
from  low  to  high  probability  of  resolution  is  usually  characterized  by  an 
increase  in  the'  false  alarm  rate.  Since  an  inordinately  large  number  of 
spurious  direction  estimates  would  be  undesirable,  the  average  number  of 
false  alarms  (per  trial)  for  each  of  the  algorithms  tested  has  been  plotted. 

Above  threshold,  the  precision  of  the  direction  estimates  is  the  primary 
measure  of  DF  performance.  In  the  statistical  literature,  an  estimate  is 
said  to  be  efficient  if  its  variance  agrees  with  the  Cramer-Rao  bound.  The 
example  in  Fig.  4.2  indicates  that  the  DF  algorithms  are  asymptotically 
unbiased  (consistent)  and  efficient  as  the  SNR  increases.  In  fact,  most  of 
the  algorithms  become  unbiased  within  a  few  dB  of  the  resolution  threshold. 
Thus,  the  remaining  issue  is  how  quickly  the  direction  error  approaches  the 
Cramer-Rao  bound.  The  simulation  data  in  Appendix  F  address  these  Issues 
over  the  range  of  parameter  values  specified  in  Table  4.1. 


43 


The  elgorithms  tested  (see  Seetlon  III.A)  were  AAR,  MEM,  MLM,  MUSIC,  and 
TNA.  Data  was  obtained  for  both  spectral  and  root  versions  of  these 
algorithms.  Since  the  root  version  of  TNA  Is  the  same  as  the  root  version  of 
AAR,  only  the  results  for  the  latter  are  presented.  In  several  Instances, 
the  false  alarm  data  for  some  of  the  algorithms  (especially  MUSIC)  nay  appear 
to  be  missing.  In  these  case,  no  false  alarms  were  observed. 

Occasionally,  the  standard  deviation  curve  for  an  algorithm  drops  below 
the  C-R  bound.  Strictly  speaking,  the  C-R  bound  plotted  In  the  performance 
comparisons  Is  valid  only  for  unbiased  and  unerpurgated  estimates.  The 
latter  reduirament  means  that.  In  every  Instance,  direction  estimates  for 
both  emitters  must  be  generated  in  order  to  Interpret  the  bound  rigorously. 
As  stated  previously,  the  DF  statistics  presented  In  this  report  are 
conditioned  on  the  emitters  being  resolved.  Near  (or  below)  threshold, 
practical  algorithms  often  tall  to  find  both  emitters  and  are  usually 
biased.  Consequently,  the  C-R  bounds  are  not  strictly  applicable  eacept  at 
high  SNR.  Moreover,  because  of  the  conditioning,  the  statistical 
significance  of  the  average  performance  decreases  as  the  probability  of 

resolution  drops* 


Be  Adaptive  Nulling  to  Support  Signal  Copy 
1.  Overview 

The  primary  function  of  adaptive  nulling  of  interference  is  to  improve 
the  signal-to-noise  interference  ratio  (SIR)  of  the  emitter  selected  for 
copy.  For  the  case  of  multiple  emitters  separated  by  less  than  one  beamwidth 
in  azimuth  angle,  the  adaptive  nulling  system  must  be  able  to  place  deep 
nulls,  within  the  "main  beam",  on  the  interferers  while  maintaining 
sufficient  gain  on  the  desired  signal  to  permit  classification  and/or 
monitoring.  In  order  to  properly  support  signal  copy,  the  nulling  system 
must  increase  the  (presumably  negative)  SIR  in  each  of  the  receiver  channels 
to  an  output  SIR  of  at  least  (say)  +10  dB  by  a  linear  weighting  of  the 
receiver  outputs. 
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As  discussed  in  Section  III. 5,  the  two  initial  candidate  nulling  methods 
under  investigation  base  their  nulling  weights  solely  upon  the  estimated 
directions  and  powers  of  all  the  emitters  present.  Therefore,  either  of 
these  nulling  techniques  is  potentially  compatible  with  any  of  the  techniques 
for  super-resolutiqn  direction-finding.  Since  there  are  several  performance 
measures  which  are  appropriate  to  the  evaluation  of  direction-finding 
algorithms  and  since  the  relative  performance  ranking  of  the  various 
algorithms  depended  upon  the  specific  performance  measure  selected,  it  was 
not  obvious  at  the  outset  which  combinations  of  DF  methods  vs.  nulling 
algorithms  would  be  most  successful  on  the  average.  Thus,  an  initial 
standardized  experiment  was  formulated  to  investigate  the  compatibility 

issue. 

The  standard  experiment  which  was  used  for  the  initial  review  of  nulling 
alternatives  was  based  upon  the  same  10-element,  half-wavelength  spaced  array 
which  was  used  for  the  initial  direction-finding  algorithm  assessment  (see 
Fig.  4.3  for  a  conceptual  block  diagram)..  Two  emitters  were  modelled  with 
the  desired  signal  -10'  dB  in  power  relative  to  the  interferer.  For  each  of 
the  candidate  pairings  of  the  DF  and  nulling  algorithms,  two  basic  parametric 
variations  were  explored:  (1)  the  effects  of  sample-size  and  (2)  the  effects 
of  residual  calibration  errors.  For  all  experiments,  the  criterion  used  to 
evaluate  the  relative  performance  was  the  same:  "what  array  signal-to-noise 
ratio  (SNR)  is  required  to  achieve  +10  dB  output  SIR  as  a  function  of  emitter 
separation?"  As  before,  array  SNR  is  measured  relative  to  the  weaker,  or  in 
this  case  desired,  signal. 

2.  Description  of  Simulation  Output 
A  typical  summary  plot,  illustrating  the  comparison  of  several 
algorithms  for  a  particular  setting  of  the  experimental  parameters,  is  shown 
in  Fig.  4.4^.  The  shaded  region  in  the  corner  of  the  plot  indicates  those 
combinations  of  array  SNR  and  signal  separation  for  which  the  indicated 
nulling  technique  (in  this  case  the  Model  Covariance  Method)  would  be  unable 
to  provide  +10  dB  output  SIR,  even  if  given  perfect  knowledge  of  all  emitter 

^The  ARM,  or  Auto  Regressive  Root  Method,  algorithm  shown  in  Fig.  4.4  is 
from  [2].  Its  performance  is  identical  with  the  root  variant  of  MEM  when 
forward— backward  averaging  of  the  covariance  matrix  is  employed. 
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Fig,  4.3.  Processing  for  simulating  adaptive  copy. 
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Fig.  4.4.  Typical  output  from  adaptive  copy  experiment. 
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directions-of-arrlval  and  received  powers.  Thus,  for  example,  at  0.1 
beamwidth  emitter  separation,  at  least  24  dB  array  SNR  of  the  target  emitter 
is  required  to  support  copy,  according  to  this  criterion. 

The  curves  in  Fig.  4.4  correspond  to  different  super-resolution 
algorithms  In  each  case,  the  DF  outputs  were  used  by  the  nulling  system  to 
set  up  the  weights  for  copy.  The  interpretation  of  the  results  is  the  same 
as  for  the  "ideal  direction-finding”  bound:  below  and  to  the  left  of  each 
curve  it  is  not  possible,  on  the  average,  for  the  indicated  DF  method  to 
provide  +10  dB  output  SIR  when  combined  with  the  Model  Covariance  Nulling 

method. 

Using  Fig.  4.4  as  a  reference  then,  we  conclude  that  about  40  dB  array 
SNR  is  required  in  order  for  MUSIC  (with  20  snapshots  of  array  data)  to 
provide  sufficiently  accurate  angle  data  on  emitters  to  permit  copy  at  0.1 
beamwidth  emitter  separation,  as  compared  with  24  dB  for  ideal  direction 
finding,  as  described  above. 

3.  Details  of  the  Experiment 

In  order  to  evaluate  the  proposed  performance  criterion,  separate  Monte 
Carlo  simulation  experiments  were  performed  for  45  pairs  of  emitter 
separations  (0.05,  0.1,  0.2,  0.4,  0.8  beamwidths)  and  array  SNR's  (10,  15, 
20,  25,  30,  35,  40,  45,  50  dB).  The  theoretically  evaluated  output  power  in 
the  true  signal  "direction,"  including  residual  calibration  errors,  if  any, 
were  tabulated.  The  desired  locus  of  +10  dB  output  SIR  was  obtained  from  a 
contour  plotting  program,  based  upon  the  rectangular  array  of  output  power 
data. 

For  each  setting  of  emitter  separation  and  array  SNR,  a  composite 
DF/copy  experiment  was  performed  as  illustrated  in  the  block  diagram  of  Fig. 
4.3.  This  diagram  depicts  the  flow  of  data  during  each  of  the  100  Monte 
Carlo  trials  over  which  the  output  SIR  was  accumulated. 

Given  the  positions  and  powers  of  the  interferer  and  signal,  the  ideal 
signal  covariance  (uncorrupted  by  calibration  errors)  is  first  computed  as 
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R  =  P  V(9  )Ae  )  +  P,  V(0^)V®(9  ) 

S  S  S  S  Li-  ■»- 


Next,  using  s  random  error  vector,  Vg,  having  Gaussian  phase  and  log-normal 
amplitude  with  the  assumed  standard  deviations,  the  steering  vectors  are 
corrupted  to  reflect  residual  calibration  effects,  such  as  unmodelled 
near-field  multipath,  etc.  This  transformation  is  obtained  by  premultiplying 
the  direction  vectors  with  a  diagonal  matrix  of  the  form 


where  aj,  <^i  are  the  n  pairs  of  amplitude  and  phase  errors  as  selected 
above.  Since  the  same  calibration  errors  are  assumed  for  both  signal  and 
interferer,  the  corrupted  signal  covariance  is  directly  obtained  by 
premultiplying  by  Re  and  post-multiplying  by  Thermal  noise  of  the 

required  intensity  is  next  added  to  the  corrupted  signal  covariance  to  yield 
the  true  covariance  of  the  signals  in  space,  as  observed  through  the  system 
with  residual  calibration  errors. 

Since  the  signal  and  the  thermal  noise  are  assumed  gaussian  conditioned 
on  the  steering  vector  errors,  the  Wishart  density  is  appropriate  for 
generating  samples,  given  the  assumed  number  of  array  snapshots.  Note  that 
this  assumes  that  the  measurement  errors,  although  unknown,  are  constant 
during  the  sampling  process.  The  sampled  covariance  matrix  is  then 
introduced  to  the  chosen  direction-finding  algorithm,  which  yields 
angles-of-arrival  and  received  powers  for  all  of  the  emitters  detected. 
Because  the  desired  signal  was  always  given  a  direction  cosine  algebraically 
lower  than  that  of  the  interferer,  the  lowest  detected  direction  cosine  from 
the  direction-finding  algorithm  was  used  by  the  copy  algorithm  as  the 
estimated  direction  cosine  of  the  signal. 
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Note  that  neither  the  DF  algorithm  nor  the  nulling  algorithm  is  given 
knowledge  of  the  corrupted  steering  vectors,  so  that  both  use  ideal, 
linear-phase,  constant-amplitude  steering  vectors.  In  order  to  determine  the 
best  nulling  weights  that  can  be  obtained  with  the  modelling  errors  present, 
therefore,  the  nulling  algorithm  is  given  the  correct  signal  DOA's  and 
powers,  but  again  is  not  told  of  the  direction  vector  errors.  The  "ideal 
direction-finding”  curve  in  the  Monte  Carlo  results  represent  a  (generally 
unachievable)  performance  limit  for  nulling  methods  which  model  the  received 
signals  with  erroneous  steering  vectors.  Only  in  the  limits  of  zero 
modelling  errors  and  infinite  array  snapshots  is  this  ideal  performance  bound 
achievable. 

The  assessment  of  the  output  SIR  is  done  with  full  knowledge  of  the 
corrupted  signal  steering  vector,  V^Og),  .and  the  theoretical 
Interference  plus  noise  covariance  matrix,  0c»  with  the  corruptions 
present.  Thus,  if  the  candidate  nulling  algorithm  yields  the  weight  vector, 
W,  then  the  output  SIR  in  the  true  signal  direction  is  computed  as  follows: 


SIR 


w^o  w 

c 


2 


The  above  quantity  is  linearly  averaged  for  the  100  Monte  Carlo  trials 
to  obtain  the  performance  statistics  for  the  overall  DF/copy  experiment. 
This  measure  may  be  Interpreted  as  the  average  performance  achievable  by  the 
calculated  nulling  vector  when  it  is  used  with  array  data  that  are 
statistically  identical  to  those  used  to  determine  the  nulling  weights. 


50 


V. 


SUTH-IARY  OF  MAJOR  FINDINGS 


A.  Direction-Finding 

One  of  the  goals  of  the  DF  investigation  was  to  determine  what  combina¬ 
tion  of  algorithms  and  system  parameters  are  needed  to  achieve  superresolu¬ 
tion,  somewhat  arbitrarily  defined  here  as  the  ability  to  resolve  emitters 
separated  by  only  0.1  beamwldth.  In  theory,  superresolution  is  possible  but, 
in  practice,  may  require  extremely  high  signal-to-noise  ratios  and/or  large 
numbers  of  looks  (snapshots).  The  simulation  results  presented  in  Appendix  F 
serve  to  quantify  these  requirements  for  a  number  of  interesting  algorithms. 

For  example,  an  examination  of  the  data  in  Fig.  5.1  indicates  that 
superresolution  is  difficult  to  achieve  with  spectral-type  algorithms  given 
only  a  modest  amount  of  data  (l.e.,  10  looks).  However,  with  more  data,  the 
situation  is  not  nearly  so  bleak.  Given  100  shapshots,  the  performance  data 
in  Fig.  5.2  show  considerable  improvement  in  most  of  the  algorithms,  par¬ 
ticularly  MUSIC.  Increasing  the  number  of  looks  by  yet  another  order  of  mag 
nitude,  the  results  in  Fig.  5.3  clearly  indicate  that  MUSIC  is  asymptotically 
much  more  sensitive  than  any  of  the  other  algorithms  tested. 

The  general  trend  of  the  data  suggests  that  the  error  in  the  MUSIC  di¬ 
rection  estimate  approaches  the  Cramer-Rao  bound  as  either  the  SNR  or  the 
number  of  looks  becomes  sufficiently  large.  In  other  words,  MUSIC  appears  to 
be  asymptotically  efficient.  However,  the  relatively  poor  sensitivity  of 
MUSIC  in  the  data-llmited  (10  snapshots)  case  is  somewhat  disappointing. 

In  Fig.  5.4,  the  10  look  experiment  has  been  repeated  with  the  emitter 
separation  increased  by  a  factor  of  two  (i.e.,  0.2  beamwidth  separation).  We 
again  see  that  spectral  MUSIC  is  noticeably  less  sensitive  than  some  of  the 
other  algorithms  (esp.,  MEM  and  AAR)  in  terms  of  its  ability  to  detect  the 
presence  of  both  emitters  at  low  SNR.  This  trend  continues  to  exist  even  at 
0.4  beamwldth  separation  (see  Fig.  5.5).  However,  in  spite  of  its  relatively 
poor  sensitivity,  MUSIC  is  generally  superior  in  terms  of  producing  more 
accurate  estimates  than  any  of  the  other  spectral  algorithms.  Of  course,  the 
SNR  must  be  sufficiently  large  to  allow  MUSIC  to  resolve  the  two  emitters. 
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Spectral  algorithms  at  0.1  beamwtdths  with  10  looks. 
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Fig*  5.2*  Spectral  algorithms  at  0.1  beamwidths  with  100  looks* 
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Spectral  algorithms  at  0,2  beamwidths  with  10 
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Spectral  algorithms  at  0,4  beamwldths  with  10  looks 


The  sensitivity  of  MUSIC  and,  in  fact,  all  of  the  algorithms  tested,  is 
improved  considerably  by  adopting  the  root-finding  approach  discussed  in  Sec¬ 
tion  III. A.  Even  when  limited  to  only  10  looks,  all  of  the  root-type  algo¬ 
rithms  are  capable  of  resolving  two  emitters  separated  by  0.1  beamwidth  at 
reasonable  signal-to-noise  ratios  (see  Fig.  5.6).  Although  MUSIC  again  ap¬ 
pears  to  be  less  sensitive  than  its  competitors,  this  slight  disadvantage  is 
relatively  insignificant  in  view  of  the  bias  results.  As  the  number  of 
available  snapshots  is  increased,  the  difference  among  the  root  versions  of 
the  DF  algorithms  become  even  less  important.  At  100  looks,  the  sensitivity 
of  all  the  algorithms  is  comparable  (see  Fig.' 5.7).  However,  MUSIC  retains  a 
discernible  advantage  in  terms  of  accuracy.  Also  significant  is  the  smaller 
bias  and  lower  false  alarm  rate  exhibited  by  MUSIC. 


B.  Adaptive  Copy 

Since  it  was  not  known  in  advance  what  the  tradeoffs  were  in  combining 
the  proposed  Model  Covariance  Method  or  the  Projection  Nulling  Method  with 
the  various  direction-finding  techniques,  all  combinations  were  tried.  For 
the  adaptive  copy  experiement  defined  in  Section  IV.B,  the  copy  algorithm 
assumed  that  the  lowest  detected  direction  cosine  was  that  of  the  desired 
signal.  Thus,  the  assumed  direction-of-arrival  could  be  far  removed  from  the 
actual  direction-of-arrival  whenever  the  direction-finding  algorithm  failed 
to  resolve  the  two  emitters.  For  this  reason,  one  might  expect  the 
direction-finding  algorithms  having  the  best  probability  of  resolution  to 
give  the  best  copy  performance.  This  was  found  to  generally  be  the  case  for 
both  spectral  and  rooting  algorithms,  as  can  be  seen  by  comparing  Fig.  5.8 
to  Fig.  F-19,  Fig.  5.9  to  Fig.  F-21,  and  Fig.  5.10  to  Fig.  F-20.  For  the 
rooting  algorithms,  this  observation  tends  to  hold  only  for  those  values  of 
array  SNR  above  that  value  where  the  roots  cross  in  angle  (a  dip  in  the 
resolution  probability  often  occurs  at  that  SNR,  as  can  be  seen  in  Fig. 

F.20). 

By  comparing  Figure  5.10  to  Figure  5.8  and  Figure  5.11  to  Figure 
5.9,  we  see  that  the  rooting  variants  of  the  direction-finding 
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algorithms  again  offer  significantly  superior  performance.  This  again 
implies  that  exploitation  of  the  array  structure  is  desired  in  order  to 
achieve  the  best  results. 

The  above  results  indicate  that  the  modeling  techniques  which  steer 
nulls  to  perform  copy  offer  nearly  optimal  performance  in  the  case  that  the 
array  is  perfectly  calibrated  for  plane  wave  direction  finding.  ■  On  the  other 
hand.  If  small  gain  and  phase  errors  are  introduced  randomly  Into  each 
receiver  channel.  Independent  of  direction-of-arrival,  there  results  a 
significant  loss  of  superresolution  copy  performance,  as  seen  in  Fig.  5.12. 
This  result  implies  that  null  steering  is  an  unacceptable  approach  for 
obtaining  superresolution  copy.  Techniques  which  exploit  the  signal  space 
decomposition  employed  in  SVD-type  direction-finding  algorithms  to  overcome 
this  deficiency  are  being  studied  as  alternatives  and  will  be  reported  in  a 
future  phase  of  this  study. 

C.  Conclusions 

Noise  cancellation  by  singular-valued  decomposition  is  a  necessary 
prerequisite  to  superresolution  (e.g.,  less  than  1/10  beamwldth). 
Exploitation  of  regular  array  structures  (e.g.,  linear  arrays)  is  necessary 
to  optimize  resolution  sensitivity  and  accuracy  of  estimates  with  small 
amounts  of  data.  ^rray  errors  due  to,  for  example,  calibration  residuals, 
cause  significant  degradation  of  superresolution  performance  and  copy. 
Methods  based  upon  open-loop  null  steering  cannot  provide  superresolution 
performance  based  upon  array  information  which  is  only  slightly  in  error 
(e.g.,  5“  phase  and/or  0.5  dB  amplitude  errors  per  element). 

Because  of  the  above  conclusions  for  the  initial  phase  of  the  study , 
focus  for  subsequent  phases  has  been  directed  toward  SVD-type  algorithms, 
both  for  direction-finding  and  copy,  and  robust  extensions  to  desensitize 
performance  to  array  errors. 
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EMITTER  SEPARATION  (Beamwidths) 


Fig.  8.  Adaptive  listening  with  spectral  superresolution 
algorithms  and  20  array  snapshots  (contours  of  10  dB  output 
SIB.  with  -10  dB  input  SIR) . 
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EMITTER  SEPARATION  (Beamwidths) 


Fig.  5.9.  Adaptive  listening  with  spectral  superresolution 
algorithms  and  100  array  snapshots  (contours  of  10  dB  output 
SIR  with  -10  dB  input  SIR) . 
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Fig.  5.10.  Adaptive  listening  with  rooting  algorithms  and 
20  array  snapshots  (contours  of  10  dB  output  SIR  with  —10 
dB  input  SIR) . 
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EMITTER  SEPARATION  (Beamwidths) 


Fig.  5.11.  Adaptive  listening  with  rooting  algorithms  and 
100  array  snapshots  (contours  of  10  dB  output  SIR  with  —10 
dB  input  SIR) . 


EMITTER  SEPARATION  (Beamwidths) 


Fig.  5.12.  Effects  of  small  calibration  errors  on  adaptive 
listening  performance  (contours  of  10  dB  output  SIR  with 
-10  dB  input  SIR) . 
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Appendix  A 

CRAMER-RAO  DIRECTION-FINDING  BOUNDS  FOR  GAUSSIAN  SIGNALS 
1,  GENERAL  DISCUSSION 

Detailed  treatments  of  the  multiparameter  Cramer-Rao  bound  are  readily 
available  in  the  engineering  literature  (e.g.,  see  [30]).  For  our  purposes, 
the  fundamental  theoretical  result  may  be  stated  as  follows: 

Consider  a  probability  density  function  p(rj  )  that  governs  a  probabil¬ 
istic  mapping  from  a  parameter  space  into  an  observation  space  r  .  The 

basic  problem  is  to  estimate  the  (vector)  parameter  0  from  the  (vector)  ob¬ 
servation  r.  The  Cramer-Rao  bound .  asserts  that  the  covariance  matrix  of  any 
unbiased  estimate  of  9  must  satisfy 

Cov(0)>f’^ 

where  F  is  the  Fisher  information  matrix  with  elements 
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mn 


=  E  { 


8  In  p 
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8  In  p 
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=  -E  { 


8^1n  p  1 

80  89 
m  n 


(A.l) 


Of  course,  the  expected  value  operation  is  defined  in  terms  of  the  given 
probability  density  function  (pdf),  i.e.. 


E  {x}  =  /  x(r)  p(rj  9)  dr 

Situations  frequently  arise  where  one  is  primarily  interested  in  a  subset  of 
the  unknown  parameters.  The  remaining  "nuisance”  parameters  are  only  impor¬ 
tant  to  the  extent  that  they  adversely  affect  estimates  of  the  "desired" 
parameters.  By  simply  partitioning  9,  i.e., 
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0  =  [  a  j  3  ]  , 

the  Fisher  matrix  may  be  written  as 

^aa  ^ae 

If  the  Fisher  matrix  has  an  inverse,  it  may  be  written  as 
t  "'aa- 

*  •  (A.2) 

^  SB  ga  aa  aB 

The  remaining  (off-diagonal)  terms  can  be  easily  obtained  but  are  not  needed 
here. 

If  the  nuisance  parameters  (e.g. ,  3)  were  known,  the  Cramer-Rao  bound 

for  an  unbiased  estimate  of  the  desired  parameters  (l.e,  a)  would  be 

Cov  (a)  >  F 

—  act 

This  inequality  is  always  valid  but  the  bound  obtained  from  (A.2)  is  tighter, 
i  •  e . , 


Cov  (i  > 

where  the  "reduced’*  information  matrix  is 

^aa  ^aa  ^aB^SS^Ba 
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In  the  following  sections,  it  will  be  convenient  to  use  the  notation 


r  =  (F  ) 

■  a  b  ab  mn 

m  n 


to  refer  to  an  arbitrary  element  of  a  Fisher  submatrix. 

2.  STATISTICAL  INFERENCE  BASED  ON  COMPLEX  GAUSSIAN  OBSERVATIONS 
Consider  a  complex  Gaussian  vector  r  with  zero  mean,  i.e., 

E  {r}  *  0 

and  circular  components,  i.e, 

E  {rr*^}  =  0 

and  covariance  R,  i#ee, 

E  {rr®}  =  R 

Note  that  a  superscript  T  indicates  the  usual  (real)  transpose  operation, 
whereas  an  H  represents  the  complex  conjugate  (Hermitian)  transpose  opera¬ 
tion.  The  unknown  parameters  9  are  imbedded  in  the  covariance,  i.e., 

R  =  R(9)  . 

Assuming  the  covariance  matrix  R  is  positive  definite,  the  complex  multivari¬ 
ate  Gaussian  pdf  [31]  is  given  by 

p(r|  9)  =  .  I'--  exp  {  -r\  ^r  }  , 
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where  L  is  the  dimension  of  r  and  jRj  < 
t reducing  the  sample  covariance  matrix 


denotes  the  determinant  of  R.  By  in- 


^  =  4  I  r(k)r®(k) 
•  ^  k=l 


the  logarithm  of  the  pdf  for  K  statistically  independent  observations  can  be 
written  as 

In  p({r(k)|k=l . k}|  e)  =  -K  r}  +  In  |e|  +  L  In  i,) 

where  Tr  {  }  is  the  standard  trace  operator.  In  arriving  at  this  result,  use 
has  been  made  of  the  Identity 

Tr  {ab}  =  Tr  {BA}  . 

In  the  next  section,  we  will  also  have  occasion  to  use 


Tr  {a”}  =  Tr*fA}  , 

Where  a  superscript  asterisk  (*)  denotes  the  complex  conjugate  operation. 

As  indicated  by  (A.l),  an  arbitrary  element  of  the  Fisher  information 
matrix  can  be  determined  by  calculating 


p  ,,  _E  (  } 

^0  6  I  39  30  J 

m  n  m  n 


Invoking  the  differential  rules  [32] 


d  R  ^  =  -r“^  dR  r“^ 


(A.5) 
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and  the  less  familiar 


d  In  |r|  =  Tr  {  r”^  dR  } 


we  first  obtain 


3  In  p 


K  Tr  {  r'^  (R”^  R  -  I)  } 

m 


Since  the  sample  covariance  matrix  is  always  an  unbiased  estimate  of  the 
actual  CO vari ance ,  i . e . , 


E  {r}  =  R 


it  follows  that 


E  { 


3^1n  p 

39  30 
m  n 


}  =  K  Tr  {  R 


-1  3R  3R 


-1 


—  R  \ 
36  30  ^ 

m  n 


^  ^  f  3R  3R~^  1 

-  K  Tr  {  je  39  J 

m  n 

Thus,  the  number  of  observations  enters  the  result  only  as  a  multiplicative 
constant.  Since  no  loss  of  generality  is  incurred,  the  simplifying  assump¬ 
tion  R=1  is  imposed  at  this  point.  Applying  (A. 5)  once  more  and  ignoring  K, 
we  finally  obtain  the  desired  result  in  a  suitable  form  for  subsequent 
calculations • 


F 


9  9 
m  n 


=  Tr  { 


3R  -1  3R  „-l  1 

36  ^  30  ^ 

m  n 


(A. 6) 
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3.  APPLICATION  TO  DIRECTION-FINDING 

For  the  direction-finding  (DF)  problem  posed  in  Section  II. A  of  this  re¬ 
port,  the  covariance  matrix  of  the  observed  vector  (snapshot)  r  can  be  writ¬ 
ten  as 

R  =  APA^  +  N 

where  the  unknown  directions 

T 

a  =  [  “l  «2  •  •  •  “j  1 

enter  the  problem  through  the  direction  matrix 
A  =  [  aCoj)  I  ...  I  aCuj)  ] 


The  only  requirement  imposed  on  the  (vector)  array  gain  a(a)  is  that  it  pos 
sess  a  derivative  in  the  usual  sense,  i.e., 

•  da 
a  = 

da 

is  assumed  to  exist  at  every  direction  a.  The  matrix  constructed  from  the 
derivatives  of  the  array  gain  at  each  of  the  unknown  directions  is  written  as 


A  =  [  a(a^)  I  ...  I  1 

The  "signal-in-space”  covariance  matrix  P  may  be  interpreted  in  terms  of  the 
pairwise  correlations  that  would  exist  among  the  J  signals  at  some  convenient 
reference  point  on  (or  near)  the  array.  Unfortunately,  the  general  problem 
where  P  is  totally  unknown  soon  becomes  overwhelming  as  the  number  of  signal 
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directions  increases.  Thus,  we  first  consider  the  ideal  case  of  known  ?.  We 
then  treat  the  important  case  where  P  is  an  unknown  diagonal  matrix.  Since 
the  signals  transmitted  by  two  emitters  are  normally  uncorrelated,  the  latter 
assumption  is  quite  reasonable  for  a  large  class  of  DF  problems. 

The  noise  covariance  matrix  N  is  assumed  to  be  non-singular  and  known. 
Under  these  conditions,  the  general  problem  is  easily  converted  via  a  linear 
transformation  to  the  special  case  where  N  is  the  identity  matrix,  i.e., 

N  — >  I 

Under  this  transformation,  the  direction  matrix  becomes 

A  — >  V  =  a 

and  similarly, 

A  — >  V  =  A  . 

The  covariance  matrix  for  the  transformed  problem  is 


R  =  VPV®  +  I 

Our  results  are  somewhat  easier  to  Interpret  if  the  Inverse  covariance  matrix 
is  expressed  in  the  form 


=  I  -  vov^ 


(A.7) 


A  direct  calculation  shows  that  0  is  the  (Hermitlan)  solution  to 


P  -  0  =  PWO  =  OWP 


(A.8) 


where  we  have  introduced  the  Gramian  matrix 


tf  =  v«v  . 

The  uniqueness  of  the  representation  in  (A. 7)  follows  from  the  fact  that 
I  +  PW  is  a  non-singular  matrix.  The  proof  is  by  contradiction,  i.e.,  sup 
pose  there  exists  a  non-zero  x  such  that 

(I  +  PW)x  =  0 

Preinultiplylng  this  equation  by  V,  it  follows  that 
(I  +  VPV®)Vx  =  0 

Since  R  is  obviously  non-singular,  it  follows  that  Vx  =  0.  In  turn,  Wx  =  0 
and  hence  the  contradiction  (x  =  0). 

Thus,  0  may  be  written  explicitly  as 

0  =  (I  +  PW)”^P  =  P(I  +  WP)“^ 

If  P  is  positive  definite,  we  also  have 

0  =  (P"^  +  W)“^ 

In  the  ensuing  development,  we  make  extensive  use  of  the  notational  device 
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for  the  partial  derivative  of  the  direction  matrix  V  with  respect  to  an  un¬ 
known  direction.  This  derivative  can  be  extracted  from  the  matrix  of  deriva¬ 
tives  defined  previously,  i.e., 

V.  =  V  e.eT  ,  (A. 11) 

.1  3  3 

where  the  unit  vector  ej  is  the  jth  column  of  the  JxJ  identity  matrix. 

Using  (A. 10),  we  first  write  the  partial  derivative  of  the  covariance 
matrix  with  respect  to  the  jth  direction  as 

-P-  =  V.PV®  +  VPV^ 

Substituting  this  expression  in  (Ae6)  and  expanding  terms,  we  get 

F  =  Tr  {  V  PV^  r”^  V  PV*^  r“^  }  +  Tr  {  V  PV®  R  ^  VPV®  R  ^  } 
a  a  m  n.  u 

m  n 

+  Tr  {  VP^  R~^  V  PV®  r"^  }  +  Tr  {  VP^  R  ^  VP^  R  •  (A.  12) 

Since  R  and  P  are  Hermitlan,  it  follows  from  the  identities  (A. 3)  and  (A. 4) 
that  the  last  term  on  the  right-hand  side  of  (A.  12)  is  the  complex  conjugate 
of  the  first  term.  Similarly,  the  second  and  third  terms  form  a  conjugate 
pair.  Further  simplifications  result  from  substituting  (A.  11)  and  arranging 
terms,  with  the  help  of  (A. 3),  to  obtain  scalar  expressions  within  the 
braces.  In  this  manner,  we  eventually  arrive  at 

F  =  2  Re  {  e^Pv\“^Ve  e^P  V®R~^Ve  +  e^Pv\~^VPe  e^v\  ^Ve  }  ,  (A.  13) 

Qtct  min  .nn  mni  n 

m  n 

where  Re  {x}  denotes  the  real  part  of  x. 
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At  this  point,  it  is  convenient  to  introduce  the  Hadamard  product  AQB 
of  two  matrices,  defined  by 

(AOB)  =  A 

'^mn  mn  ran 

for  any  two  matrices  with  the  same  dimensions.  Unlike  the  usual  matrix 
product,  the  Hadamard  product  is  commutative,  i.e.,  A[]B  =  B[]A.  Since 

(a'^DB)_  =  A  B 
^  tnn  tun  nm 

the  Fisher  (sub)matrix  with  elements  specified  by  (A. 13)  can  be  written  as 

F  =  2  Re  {(Pv“r‘^V)^0(PV°r‘^V)  +  (PV®r‘^VP)^D  (^R"^V)  }  (A.  14) 

act 

The  matrix  expressions  that  appear  in  this  result  can  be  simplified 
considerably  by  using  (A.7),  (A.8),  and  (A. 9),  e.g., 

PV\"^VP  =  PV^d  -  V0V®)VP 
=  (P  -  PWO)WP 

=  0\TP 


Similarly, 

pv®r“^v  =  PV^d  -  vov’*)  V 
=  (P  -  PWO) 
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where,  in  the  last  step,  we  have  introduced 


w  =  v®v 

despite  the  apparent  abuse  of  notation.  Consequently,  we  may  also  write 

v\“^V  =  v”(I  -  VQV”)  V 
=  -  H^QW- 

Substituting  these  identities  in  (A. 14),  the  final  form  of  the  Fisher  matrix 
for  the  unknown  directions  a  emerges  as 

F  =  2  Re  {(P-O)^O(^  -  W^OW)  +  (OW)'^D(QW)  } 
act 

In  order  to  proceed  without  undue  difficulty,  we  now  assume  that  the 
signal-in-space  covariance  P  is  a  positive  (definite)  diagonal  matrix  and, 
accordingly,  introduce  the  nuisance  parameters 

-  1" 

The  partial  derivatives  of  the  covariance  matrix  R  with  respect  to  the 
nuisance  parameters  are  given  by 


It  follows  easily  from  (A. 6)  that 


* 


T  n  — 1  T  H  -1 

P  Ve  e  V  R  P  Ve  e  V  R 
mm  m  m  nn  n  n 


=  P  e'^V®R~^Ve 
Tinn  m  n. 


P  e'^V^R  ^Ve 
nn  n  m 


Since 


(DA) 


mn 


=  D  A 
'  mm  mn 


for  any  diagonal  matrix  D,  we  may  write  the  submatrix  for  the  nuisance 
parameters  as 


=  (PV®R  ^V)  □  (PV^'^V)*^ 

PP 


From  the  development  in.  (A.  15),  we  deduce  the  identity 


PV^R~^V  =  CW 


(A. 17) 


which  yields 


Fgg  =  (OW)  □  (OW)’’^ 


The  calculation  of  the  coupling  matrix  proceeds  along  similar  lines,  e.g. , 


S  a 
m  n 


=  Tr  {  R 


-1  3R  ^-1  3R  1 

36  ^  3a  * 

m  n 
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=  Tr  {  P  R"^Ve  (V  PV^  +  VPV®)  } 

^  mm  m  m  n  ri 


=  2  Re  {  P  e^v\"^Ve  e^Pv\"^Ve  } 

*•  mm  m  n  ti  m 

Consequently,  the  information  coupling  between  the  directions  and  the 
nuisance  parameters  is  determined  by 

=  2  Re  {  (PV®R"^V)  □  (PV^r’^V)*^  } 

pa 

Simplifying  with  (A. 16)  and  (A.  17)  yields 

=  2  Re  {  (OW)  0  (QW)^  } 
pa 

Finally,  since  a  Fisher  information  matrix  is  symmetric,  it  follows  that 
T 

F  «  =  Ft 
a6  Sa 
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Appeadix  B 

LINEAR  PREDICTION  AND  MAXIMUM  ENTROPY  SPECTRAL  ANALYSIS 


Maximum  entropy  spectral  analysis  is  closely  linked  to  one  of  the  more 
fascinating  questions  in  time  series  analysts  —  given  a  record  of  the  past, 
how  well  can  one  predict  the.  future?  Differences  of  opinion  concerning  such 
matters  are  what  horseracing  is  all  about;  however,  within  the  relatively 
orderly  realm  of  stationary  random  processes,  the  theory  of  linear  prediction 
offers  a  fairly  complete  answer.  The  discussion  in  this  appendix  covers  most 
of  the  salient  points  of  this  theory  for  complex-valued  processes. 

1.  LINEAR  PREDICTION 

Let  e(n)  denote  the  minimum  mean  square  error  that  can  be  achieved  by 
linearly  combining  the  n  most  recent  observations  of  a  zero-mean,  stationary 
complex  random  process  in  order  to  predict  its  next  value.  Intuitively,  one 
would  expect  the  prediction  error  to  decrease  as  the  number  of  available  ob¬ 
servations  increases.  Indeed,  it  is  a  well-known  fact  that  the  prediction 
error  satisfies  a  recursive  relationship  of  the  form 

e(n)  =  [  1  -  jk(n)|^l  e(n-l)  (B.l) 

where  k(n)  is  the  nth  (complex)  reflection  coefficient  of  the  process.  Since 
the  prediction,  error  is  always  non— negative,  the  magnitude  of  a  reflection 
coefficient  can  never  exceed  unity.  Unless  otherwise  explicitly  stated,  per¬ 
fectly  predictable  (i.e.,  degenerate)  processes  are  excluded  from  considera¬ 
tion  in  order  to  obtain  the  strict  inequality 

|k(n)|  <  1  . 

Makhoul  [10]  has  emphasized  the  equivalence  between  the  reflection  coeffi¬ 
cients  of  a  process  and  its  correlation  coefficients.  Unless  a  process  is 
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degenerate,  the  first  n  correlation  coefficients  uniquely  determine  the  first 
n  reflection  coefficients  and  vice  versa.  This  statement  also  holds  for  de¬ 
generate  processes  provided  e(n-l)  >  0. 

Reflection  coefficients  provide  an  extremely  convenient  way  to  generate 
prediction  filters  recursively.  Let  the  vector  w(n)  denote  the  weights  of 
the  nth  order  prediction  error  filter  based  on  the  n  most  recent  observa¬ 
tions.  As  shown  in  Fig.  B.l,  a  prediction  error  filter  is  constructed  by 
subtracting  the  output  of  the  corresponding  prediction  filter  from  its  in¬ 
put.  Given  the  filter  weights  obtained  at  the  (n-l)st  stage  of  the  recur¬ 
sion,  the  nth  order  weights  are  completely  determined  by  the  nth  reflection 
coefficient.  This  relationship  may  be  stated  as 

w(n)  =  w_j_(n-l)  -  k(n)  J  w*(n-l)  (B.3) 

where  the  subscript  "+"  denotes  a  vector  which  has  been  extended  by  appending 
a  trailing  zero,  e.g., 

[a  b]^  =  [a  b  0]  , 

and  J  is  the  appropriate  exchange  matrix  that  simply  reverses  the  order  of 
the  elements  within  a  vector.  For  example,  the  3x3  exchange  matrix  is 

0  0  1 
0  10 
1  0  0 

and 

[a  b  c]  J  =  [c  b  a] 

The  self-inverting  property  of  exchange  matrices,  i.e.. 
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Is  often  useful. 

The  prediction  error  filter  of  order  zero  has  no  memory,  and  the  recur¬ 
sion  in  (B.3)  begins  with  w(0)  =  1.  Actually,  the  leading  element  of  w(n)  is 
always  unity  —  a  consequence  of  the  direct  signal  path  in  the  block  diagram 
defining  the  prediction  error  filter  (see  Fig.  B.l).  A  more  noteworthy  ob¬ 
servation  is  that  the  last  coefficient  of  the  nth  order  prediction  filter  is 
the  nth  reflection  coefficient.  This  interesting  property  enables  one  to  re¬ 
construct  the  first  n  reflection  coefficients  from  the  nth  order  linear  pre 
diction  (LP)  filter.  The  procedure  for  accomplishing  this  task  is  essen¬ 
tially  the  same  as  the  famous  Schur  stability  test  [33].  The  key  here  is  to 
realize  that  w(n-l)  can  be  obtained  directly  from  w(n),  e.g., 

w_j_(n-l)  =  (1  -  |k(n)|^)“^  [w(n)  +  k(n)Jw*(n)] 


follows  easily  from  (B.3) 

Of  course,  the  coefficients  for  a  prediction  error  filter  may  be 
obtained  directly  by  minimizing  its  expected  output  power.  In  order  to 
pursue  this  (equivalent)  approach,  let  the  vector  r  represent  an  arbitrary 
segment  of  the  observed  process.  The  covariance  (correlation)  matrix  of  r  is 
usually  written  as 

R  =  Efrr^} 

where  a  superscript  H  denotes  the  complex  conjugate  (Hermitian)  transpose 
operation.  The  response  of  a  finite  impulse  response  (FIR)  filter  to  a  seg¬ 
ment  of  the  process  can  be  written  as 

Tr 

y  =  w  Jr  , 
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where  w  represents  the  coefficients  (impulse  response)  of  the  filter,  and  the 
superscript  T  denotes  the  usual  (real)  transpose  operation. 

The  Tnlnitnum  mean  square  prediction  error  is 

e  =  min  E{jyj^} 

T^t  *  (B.5) 

=  min  w  JRJ  w  , 

where  the  minimization  is  over  all  w  with  a  leading  coefficient  of  unity. 

When  the  observed  process  is  statistically  stationary,  R  is  a  Toeplitz 
matrix  of  (complex)  correlation  coefficients,  i.e.. 


(B.6) 

^m-n 

Since  a  covariance  matrix  is  always  Hermitian,  the  correlation  coefficients 
(lags)  must  also  satisfy 

*  (B.7) 

=-n  =  ^  n  • 

Finally,  without  suffering  any  significant  loss  of  generality,  we  may 
restrict  our  attention  to  normalized  processes  with  unit  variance,  i.e., 

Cg  =  1 

Given  (B.6)  and  (B.7),  it  is  a  relatively  trivial  exercise  to  show  that  R  has 
the  (harmonic)  property 
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(B.8) 


JRJ  =  R* 

Substituting  (B.8)  in  (B.5),  the  linear  prediction  problem  may  be  stated  as 
follows.  Find  the  filter  weights  which  achieve  the  minimum  mean  square  error 

e  =  min  w^w  , 
subject  to  the  constraint 
II  0  ...  0]  w  =  1 

This  constrained  minimization  problem  is  easily  solved  using  the  method  of 
Lagrangian  multipliers.  A  complete  solution  is  obtained  by  solving  the 
(mixed)  linear  system  of  equations 

Rw  =etl0...0]^  .  (B.9) 

Thus,  the  first  n  lags  of  a  stationary  process  determine  the  nth  order  pre¬ 
diction  filter.  The  following  argximent  establishes  the  converse. 

Using  the  method  of  bordering,  it  can  be  shown  [34]  that  (B.l)  and  (B.3) 
solve  (B.9)  recursively.  '  For  non— degenerate  processes,  the  nth  reflection 
coefficient  is 

k(n)  =  c^(n)Jw^(n-l)  /  e(n-l)  (B.IO) 


where  the  vector 

T 

c(n)  =  [1  Cj  ...  c^] 
is  constructed  from  the  first  n  lags. 
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Given  a  prediction  error  filter  w,  we  may  (now)  invoke  <B.4)  to  generate 
all  the  lower  order  prediction  error  filters.  Following  Burg  [35],  we  con¬ 
struct  a  lower  triangular  matrix  L,  column  by  column,  from  the  weights  for 
the  error  filters  of  order  n  through  0.  Since  the  first  n  reflection  coef¬ 
ficients  can  be  extracted  from  L,  e.g., 

[0  ...  0  1]  L  =  I-k(n)  ...  -k(l)  1] 

we  may  compute  the  corresponding  prediction  errors  from  (B.l),  starting  with 
e(0)  =  1.  The  prediction  errors  are  placed  along  the  main  diagonal  of  a 
( (ji agonal )  matrix  D  in  order  of  increasing  value. 

By  taking  advantage  of  the  Toeplitz  structure  of  the  correlation  matrix 
R,  the  prediction  error  filter  Eq.  (B.9)  can  be  generalized  to 


where  U  is  an  upper  triangular  matrix  (to  be  determined).  By  construction, 
the  main  diagonals  of  L  and  IT  consist  of  all  "r's.  Pre-multlplylng  (B.ll)  by 
the  Hermitian  transpose  of  L,  we  obtain 


L®RL  = 


Since  the  product  of  two  lower  (upper)  triangular  matrices  is  also  a  lower 
(upper)  triangular  matrix,  the  left-hand  side  of  (B.12)  is  both  an  upper  and 
a  lower  triangular  (Hermitian)  matrix,  i.e.,  a  real  diagonal  matrix.  It  fol- 
lows  that 

A  =  u\ 
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is  a  real  diagonal  matrix  with  "l"s  along  the  main  diagonal  (i.e.,  an  iden¬ 
tity  matrix).  Therefore,  U  is  the  inverse  of  the  Hermitian  transpose  of  L, 
and  consequently,  the  correlation  matrix  R  is  given  by 

R  =  =  OT)U^  .  (B.13) 

Since  L  and  D  were  constructed  from  the  prediction  error  filter  w,  we  may 
conclude  that  the  nth  order  LP  filter  for  a  stationary  process  uniquely  spec¬ 
ifies  the  first  n  lags  of  the  process. 

If  the  determinant  of  the  nxn  correlation  matrix  R  is  written  as  d(n), 

it  follows  from  (3.13)  that 


n 

d(n+l)  =  n  e(k) 
k=0 


or  equivalently, 

e(n)  =  d(n+l)/d(n)  .  (B.14) 

2.  MAXIMUM  ENTROPY  SPECTRAL  ANALYSIS 

The  "ratio  of  determinants”  formula  ultimately  enables  us  to  relate  the 
error  in  predicting  a  process,  given  its  entire  past,  to  the  entropy  of  the 
process.  Taking  the  limit  of  (3.14)  for  arbitrarily  large  n  and  applying  a 
well-known  theorem  from  real  analysis  [36]  yields 

e(oo)  =  iim  e(n) 

=  lim  d(n+l)/d(n) 

=  lim  2.n  d(n)l^^^ 
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Taking  the  logarithm  of  both  sides  leads  to 

In  e(«)  =  lim  (l/n)  In  d(n)  .  (B.15) 

Although  it  is  not  essential  for  this  discussion,  the  right-hand  side  of  this 
result  fan  be  expressed  in  closed  form,  e.g.,  see  [37], 


lim  (l/n)  In  d(n)  -  (l/2r) 


/ 


-ir 


In  SCoj)  do. 


(B.16) 


where 


S(a,)  =  l  +  2  ^  c  cos  no. 

n=l  n 

is  the  true  power  spectral  density  of  the  process  (i.e.,  the  Fourier  series 
of  the  correlation  coefficients).  The  absolute  integrability  of  In  S((l)  is 
known  in  the  literature  as  the  Paley-Wiener  condition.  The  right-hand  side 
of  (B.16)  is  sometimes  referred  to  as  the  entropy  (rate)  of  the  power  spec¬ 
tral  density  S(a,). 

The  entropy  of  a  random  vector  r  with  probability  density  function  p(r) 
is  defined  to  be 

h  =  E  {  -In  p(r)  } 

In  particular,  the  entropy  of  an  n-dlraensional  circular  (complex)  Gaussian 
vector  with  covariance  R  is  given  by  [see  Appendix  A,  section  2). 

h(n)  =  n  (1  +  In  ir)  +  In  |r|  .  (B.17) 

Thus,  the  average  entropy  per  observation  of  a  stationary  complex  Gaussian 
process  is 
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<h>  =  lim  (1/n)  h(n) 


=  1  +  In  IT  +  llm  (1/n)  In  d(n) 

Setting  n  =  1  in  (B.17),  we  note  that  the  entropy  of  a  circular  Gaussian 
random  variable  with  unit  variance  is 


h  =  h(l) 


Substituting  (B.15)  and  (B.19)  in  (B.18),  the  relationship  between  the 
entropy  (rate)  of  a  complex  stationary  Gaussian  process  and  its  optimum 
prediction  error  emerges  as 


<h>  =  h  +  In  e(“) 


(B.20) 


Thus,  the  greater  its  entropy  (i.e.,  disorder),  the  more  difficult  a  process 
is  to  predict. 

At  this  point,  we  are  to  consider  the  class  of  processes  that  share  the 
same  given  set  of  N  initial  lag  values.  Of  course,  every  process  in  this 
class  has  the  same  prediction  error  based  on  the  N  most  recent  observations. 
Given  additional  (older)  observations,  the  worst  conceivable  situation  is  for 
the  prediction  error  to  remain  constant.  The  implication  of  this  rather 
pessimistic  assumption  is  that  the  additional  observations  are  absolutely 
useless. 

Pessimistic  or  not,  the  situation  described  above  is  a  legitimate  model 
and  corresponds  to  choosing  the  (unknown)  reflection  coefficients  to  be 
identically  zero  for  n  >  N.  For  this  model,  we  have 


e(»)  =  e(N) 
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and,  in  light  of  (B.20),  the  average  entropy  is  clearly  maximized.  Moreover, 
(B.3)  shows  that  the  Nth  order  maximum  entropy  model  has  the  property  that 
its  optimum  prediction  filter,  based  on  the  infinite  past,  is  the  same  as  its 
Nth  order  prediction  filter. 

One  of  the  most  important  consequences  of  the  orthogonality  principle 
[38]  is  that  the  output  of  an  optimum  prediction  error  filter  is  a  white 
noise  process.  Consequently,  a  statistically  equivalent  input  process  can  be 
generated  by  driving  the  inverse  of  the  optimum  error  filter  with  white 
noise.  Referring  to  Fig.  B.l,  elementary  linear  systems  analysis  shows  that 
the  required  Inverse  can  be  constructed  using  the  simple  feedback  arrangement 
depicted  in  Fig.  B.2.  Since  the  maximum  entropy  error  filter  has  a  finite 
memory  (impulse  response),  its  inverse  is  an  all-pole  filter.  Thus,  the 
critical  assumption  behind  the  maximum  entropy  method  is  that  the  observed 

time  series  is  an  autoregressive  process. 

The  validity  of  the  autoregressive  assumption  is  an  issue  far  beyond  the 
scope  of  this  report.  Dltimately,  a  preference  for  one  model  over  another 
should  be  based  on  empirical  considerations.  Deliberations  about  such  mat¬ 
ters  could  perhaps  be  made  more  meaningful  by  the  Intelligent  application  of 
statistical  hypothesis  testing  techniques  [30]. 

The  point  to  be  made  here  is  that  the  maximum  entropy  method  (MEM)  al¬ 
ways  leads  to  legitimate  spectral  estimates.  An  MEM  spectral  estimate  takes 

the  form 


N 


1  +  w  exp(-ino.)  I' 
iH-1  " 


(B.21) 


where  the  coefficients  are  obtained  from  the  underlying  LP  filter 


w  =  [1  w, 


• 
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and  e  is  the  corresponding  prediction  error.  In  the  traditional  Blackman- 
Tukey  approach,  a  spectral  estimate  is  derived  from  N  lag  estimates  by  com¬ 
puting  the  (truncated)  Fourier  series 

N 

S(tt,)  =  1  +  2  ^  c  cos  n(jb  .  ' 

n+1 

Unfortunately,  this  estimate  is  not  generally  guaranteed  to  be  positive! 

The  difficulty  with  the  latter  approach  can  be  attributed  to  the  fact 
that  the  calculation  in  (B.22)  implicitly  assumes  that  the  (unknown)  lags  for 
n  >  N  are  zero.  The  standard  "fix”  for  this  problem  is  to  pre-multiply  the 
lags  by  a  window  function  that  is  zero  for  n  >  N  and  has  a  positive  Fourier 
series  (transform).  This  class  of  functions  has  been  studied  extensively 
[39],  and  a  window  can  usually  be  found  that  yields  a  positive  spectral  esti¬ 
mate*  However,  a  spectrum  computed  from  "windowed"  lags  is  never  entirely 
consistent  with  the  original  Information  (lag  values). 

Burg  was  disturbed  by  the  distortion  introduced  by  a  more  or  less  arbi¬ 
trarily  selected  window  function,  and  he  argued  that  the  given  set  of  corre¬ 
lation  coefficients  should  be  extended  to  produce  the  spectral  estimate.  Of 
all  possible  extensions.  Burg  preferred  the  extension  with  maximum  entropy. 
As  we  have  already  seen,  this  choice  amounts  to  truncating  the  reflection  co¬ 
efficients  rather  than  the  lags.  One  may  certainly  question  whether  maximiz¬ 
ing  entropy  is  the  "correct"  approach,  but  the  fact  remains  that  the  maximum 
entropy  method  always  yields  a  positive  spectral  estimate  (B.21)  consistent 

with  the  original  lags. 

3.  A  MOTE  ON  DEGENERATE  PROCESSES 

If  the  process  is  perfectly  predictable  (i.e.,  degenerate)  after  M  ob¬ 
servations,  the  LP  filters  are  no  longer  uniquely  specified  for  n  >  M.  Faced 
with  this  situation,  one  may  prefer  the  filter  with  minimum  norm.  The  coef¬ 
ficients  for  this  filter  are  easily  computed  by  simply  choosing  the  nth  re¬ 
flection  coefficient  to  minimize 
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jw(n)|^  = 

for  n  >  M»  The 


k  . 

mn 


(n) 


jw_|_(n-l)  -  k(n)  Jw_|_(n“l)j 
reflection  coefficient  that  minimizes 
w^(n-l)  J  w_^(n-l) 


this  expression  is 


In  this  case,  the  Schwartz  inequality  guarantees  that  the  magnitude  of  the 
reflection  coefficient  is  less  than  unity.  Filter  weights  calculated  via  the 
minimum  norm  criterion  obey  the  (now)  familiar  recursive  relationship 

|w(n)|^  =  [1  -  |^nin^“>|^l  |w(n-l)| 


91 


Appendix  C 

COMPARISON  OF  ALTERNATIVE  LIKELIHOOD  RATIO  TESTS 
FOR  ESTIMATING  THE  NUMBER  OF  SIGNALS  PRESENT 

As  was  discussed  earlier.,  the  MUSIC  algorithm  needs  an  estimate  for  the 
number  of  emitters  in  order  to  form  the  direction  finding  spectrum.  In  this 
appendix,  we  present  the  results  of  a  comparative  study  of  several  candidate 
likelihood  ratio  tests  (LRT’s)  which  have  origins  in  the  multivariate 
statistics  literature  [42-45].  Simkins  [46]  and  Schmidt  [47]  have  proposed 
the  use  of  an  LRT  for  driving  the  MUSIC  algorithm.  The  number  of  signals  is 
estimated  as  follows:  A  sequence  of  likelihood  ratios  is  formed  using  the 
eigenvalues  of  the  sample  covariance  matrix  for  the  hypotheses  of  zero 
signals,  one  signal,  etc.  Their  values  are  compared  to  a  user-selected 
confidence  level  of  the  appropriate  distributions,  and  the  first 

hypothesis  whose  likelihood  ratio  passes  the  test  is  accepted. 

1.  Candidate  LRT's  for  use  with  MUSIC 

Figure  C.l  illustrates  the  direction-finding  accuracy  achieved  by  MUSIC 
when  driven  by  various  LRT’s.  ^  As  a  benchmark,  ,  the  performance  of  MUSIC 
when  the  ntamber  of  signals  is  assumed  to  be  two  is  also  shown  (as  "MUSIC  w/o 
LRT").  The  best  performance  was  achieved  by  MUSIC  driven  by  the  Simkins  LRT 
using  the  Lawley  approximation  [42-44]^  -  it  was  able  to  achieve  the 

benchmark  standard.  All  the  LRT’s  studied  asymptotically  approach  a 
distribution  [42-46];  the  Lawley  approximation  helps  to  match  the  LRT  to  the 


iFor  Figs.  C.l  and  C.3,  the  statistics  were  accumulated  based  on  the 
condition  that  at  least  one  signal  was  detected.  For  assessment  purposes, 
when  only  one  signal  was  resolved,  the  angle  estimate 

undetected  signal  was  set  equal  to  that  of  the  first.  At  low  SNR  s  the  .lUSI 
spectrum  exhibits  one  peak  located  at  the  centroid  of  the  two  emitters.  As  a 
result,  the  RMS  error  tends  to  equal  half  the  emitter  separation  for  the  case 

of  equal“powered  emitters. 


%ie  form  for  the  extra  approximation  term  given  by  Simkins  [44]  is 
incorrect,  and  should  be  replaced  by  the  form  given  by  Lawley  [41]. 
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rms  ERROR  (Beamwidths) 


Fig.  C.l.  Comparison  of  MUSIC  accuracies  obtained  with  various  likelihood 
ratio  tests. 
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assumed  distribution  for  small  numbers  of  snapshots  [43,  46].  Without  the 
Lawley  approximation,  the  LRT  overestimates  the  number  of  signals.  As  a 
result,  "MDSIC  with  the  Simkins  LRT"  performs  poorer  than  "MUSIC  with  the 
Simkins  LRT  using  the  Lawley  approximation".  Finally,  "MUSIC  with  the 
Schmidt  LRT”  is  not  at  all  matched  well  to  the  problem.  The  Schmidt  LRT  [47] 
exhibits  poor  performance  because  the  distribution  being  used  has  the 

wrong  number  of  degrees  of  freedom  for  the  complex  valued  sample  covariance 

matrices  under  consideration. 

2.  Candidate  LRT's  for  use  with  ROOT  MUSIC 

Because  of  their  poor  performance  in  the  MUSIC  LRT  comparison,  the 
Schmidt  LRT  and  the  Simkins  LRT  (without  the  Lawley  -approximation)  .were 
automatically  ruled  out  for  use  with  ROOT  MUSIC.  A  comparison  of  ROOT  MUSIC 
driven  by  the  Simkins  LRT  with  the  Lawley  approximation  versus  ROOT  MUSIC 
assuming  two  emitters  (i.e.,  "w/o  LRT"),  as  seen  in  Figs  C.2  and  C.3, 

demonstrated  a  need  for  a  more  sensitive  LRT. 

All  of  the  previously  considered  LRT's  make  no  assumption  regarding  the 
noise  power.  For  the  problem  at  hand,  we  decided  that  it  was  reasonable  to 
assume  a  known  noise  level,  so  we  tested  an  LRT  which  takes  advantage  of  this 
knowledge  —  what  we  refer  to  as  the  "Lawley  LRT  (known  noise  power)  [43]. 
As  seen  in  Fig.  C.4,  the  Lawley  LRT  detects  the  presence  of  two  emitters 
approximately  4  dB  earlier  than  the  Simkins  LRT  with  the  Lawley 
approximation, resulting  in  improved  performance  in  the  ROOT  MUSIC  algorithm 
(Figs.  C.2  and  C.3).  Note  that  in  Fig.  C.2  it  is  apparent  that  perhaps  an 
additional  5  dB  of  sensitivity  might  be  obtainable,  but  that  its  impact 
(Fig.  C.3)  upon  estimation  performance  would  not  be  significant.  Thus,  the 
sinking  LRT  with  the  Lawley  approximation  was  used  for  ROOT  MUSIC  Monte  Carlo 
experiments  described  in  Section  IV. 

3 .  Summary 

The  Simkins  LRT  with  the  Lawley'  approximation  was  judged  adequately 
sensitive  to  drive  the  MUSIC  algorithm,  regardless  of  the  number  of  snapshots 
available.  The  Lawley  LRT  for  known  noise  power  provides  needed  additional 
sensitivity  to  drive  ROOT  MUSIC  and,  of  course,  could  also  be  used  with 
MUSIC.  There  remains  a  study  of  the  sensitivity  of  the  performance  of  the 
Lawley  LRT  to  errors  in  estimating  the  noise  power. 
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PROBABILITY  OF  DETECTION 


Fig.  C.2.  Comparison  of  ROOT-MUSIC  resolution  performance  with  various 
likelihood  ratio  tests. 
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RMS  ERROR  (Beamwidth) 


Fig,  C.3,  Comparison  of  ROOT-MUSIC  accuracies  obtained  with  various 
likelihood  ratio  tests. 
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LAWLEY  LRT 
(Known  Noise  Power) 


SIMKINS  LRT  USING 
-  LAWLEY  APPROX. 
(Unknown  Noise  Power) 


Fig.  C.4.  Detector  output  comparisons 


for  two  likekihood  ratio  tests 


Appendix  D 

COMPARISON  OF  TWO  POWER  ESTIMATION  TECHNIQUES 
VIA  COMPUTER  SIMULATION 

An  angle  estimation  algorithm  produces  a  set  of  candidate  signal 
directions  and  tests  them  for  authenticity  by  computing  a  power  estimate  for 
each  one.  To  obtain  a  set  of  direction  estimates,  a  vector  of  weights  must 
be  computed  which  minimizes  the  total  power  output  of  the  array  under  various 
constraints  which  serve  to  prevent  the  weighting  vector  from  being  the  zero 
vector.  Many  such  constraints  are  possible,  leading  to  many  different 
direction  finding  algorithms.  With  total  power  minimized,  the  gain  in  each 
of  the  true  signal  directions  must  also  be  minimized.  Therefore,  the 
locations  of  nulls  in  the  weighted  array  pattern  (or  peaks  in  the  reciprocal 
pattern)  provide  the  desired  set  of  candidate  directions.  An  example  of  a 
plot  of  the  reciprocal  array  pattern  produced  by  the  MUSIC  algorithm  is  shown 
in  Fig.  D.l  for  a  uniform  linear  array  of  five  elements  spaced  X/2  apart. 
There  are  two  signals  located  at  0  and  2/3  (sin6  coordinates) ,  and  since 
there  are  no  ambiguities  in  this  array,  the  two  peaks  accurately  reflect  the 
directions  of  the  signals.  Notice  that  these  are  the  only  two  peaks 
appearing  in  the  reciprocal  pattern# 

Unfortunately,  inany  angle  estimation  algorithms  produce  more  candidate 
directions  than  the  number  of  actual  signals  present,  with  the  true  signal 
directions  scattered  among  spurious  candidate  directions.  This  is  especially 
true  for  thinned  linear  arrays  due  to  the  existence  of  ambiguities.  An 
ambiguity  occurs  when  the  direction  vectors  associated  with  a  set  of 
directions  are  linearly  dependent.  True  signals  located  in  these  linearly 
dependent  directions  will  cause  false  signal  indications  in  other  directions 
which  are  linearly  dependent  upon  the  directions  of  the  actual  signals. 

To  more  clearly  illustrate  this  concept  of  linearly  dependent 
directions,  refer  to  Fig.  U.2.  This  shows  how  ambiguous,  linearly  dependent 
angles  are  created  when  the  center  element  of  the  aforementioned  five-element 
uniform  array  is  removed.  For  emitters  at  the  angles  of  ±41.8”  (sin  41.8  = 
2/3)  to  the  thinned  array,  the  phase  at  each  element  is  shifted  ±120®, 
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SPECTRUM  HEIGHT  (dB) 


SINE  OF  OFF-BROADSIDE  ANGLE 


Fig.  D.l.  MUSIC  spectrum  with  infinite  looks  for  a  uniform  five-element 
array . 
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13  3937-1 


EMITTER 
AT  0° 

EMITTER  I 

AT  -42° 


Fig.  D.2.  Example  of  linearly  dependent  direction  vectors  for  thinned 
five“element  array. 
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combined  with  an  emitter  at  0°,  the  total  signal  at  each  element  adds  to 
zero.  The  signal  does  not  add  to  zero  at  the  center  element,  but  since  that 
element  has  been  removed,  this  information  is  not  available.  Due  to  the 
linear  dependence  of  the  direction  vectors,  signals  located  at  0°  and  42° 
would  be  indistinguishable  from  a  pair  located  at  ±42°,  or  0°  and  -42  . 
Furthermore,  this  ambiguity  is  amplitude  independent  and  occurs  for  any  set 
of  three  angles  which  satisfy 

sine^  -  sln02  =  sin 02  -  sln0j^  =  2/3 

Figure  D.3  is  the  reciprocal  pattern  plot  from  the  four-element  thinned  array 
mentioned  above.  The  true  signals  are  still  located  at  0  and  2/3  as  in  the 
other  pattern,  but  notice  how  removing  the  center  antenna  element  causes  the 
angle  estimation  algorithm  to  show  a  third  spurious  peak  at  -2/3. 

The  purpose  of  the  power  estimation  algorithm  is  now  clear.  Hopefully, 
if  a  signal  power  estimate  is  computed  for  each  candidate  signal  direction, 
then  the  powers  in  the  spurious  directions  will  be  at  the  noise  level  or 
significantly  below  those  of  the  true  signal  directions.  A  simple  threshold 
test  can  then  identify  and  reject  the  "false  alarms"  while  retaining  the  true 
signal  directions. 

The  most  straightforward  way  to  estimate  the  power  of  each  signal  is  via 
a  direct  mathematical  solution.  First,  assume  that  the  sampled  covariance 
matrix,  S,  is  equal  to  the  true  covariance  matrix,  R,  which  is  unknown.  It 
may  be  shown  that  any  true  covariance  matrix  can  be  written  in  the  form 

R  =  a^I  +  VPV^ 

where  the  columns  of  V  are  the  direction  vectors  and  P  is  a  diagonal  matrix^ 
whose  entries  correspond  to  the  power  for  each  direction  vector  v.  Solving 
the  above  equation  for  P  with  S  substituted  for  R  yields 


the  signals  are  incoherent. 
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SPECTRUM  HEIGHT  (dB) 


Fig.  D.3.  MUSIC  spectmjm  with  infinite  looks  for  a  thinned  linear  array. 
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p  =  (V®V)"^  (S  -  oh)  V  (V^V)”^ 

When  P  is  determined  in  this  manner,  it  will  not  be  exactly  diagonal, 
since  S  =  R  was  only  an  approximation.  The  diagonal  elements  of  P  are 
nonetheless  a  very  reasonable  approximation  of  the  signal  powers,  so  the 

non— diagonal  elements  may  be  ignored. 

The  direct  solution  method  will  work  unless  (V^V)  is  singular.  This 
occurs  when  the  columns  of  V  are  dependent.  A  singular  (V^V)  is  not  a 
problem  for  uniform  linear  arrays  with  interelement  spacing  <  Xf2  because  no 
direction  vectors  are  dependent,  but  it  is  a  problem  for  thinned  arrays. 

To  avoid  this  problem  a  least-squares  solution  is  sought  in  which  the 
powers  Pi,  P2,  •  •  •»  Pn  are  selected  so  as  to  minimize  the  squared 

error  between  S  and  R.  The  squared  error  may  be  expressed 

E  =  ns  -  Rll-  (Frobenlus  norm) 

R 


I  I 

m  n 


Since  the  inverse  of  V^V  is  not  required  to  solve  this  problem,  this 
method  has  a  chance  of  exhibiting  reasonable  performance  when  signals  are 
located  ambiguously.  In  an  ambiguous  case,  the  actual  signals  should  have 
normal  power  estimates  and  the  spurious  signal  indications  should  have  power 
estimates  at  the  noise  level. 

In  order  to  compare  the  theoretical  performance  of  the  direct  solution 
and  the  least-squares  solution  in  an  ambiguous  situation,  a  computer 
simulation  was  set  up  using  the  four-element  thinned  array  as  the  receiving 
antenna  (Fig.  D.4)  and  two  emitters  of  equal  power.  One  emitter  was  held 
fixed  at  0®  or  directly  broadside  to  the  antenna.  The  other  emitter  was 
moved  step  by  step  closer  to  41.8“  (Arcsln  2/3)  which  is  an  ambiguous 
situation.  The  bearawidth  of  this  antenna  is  23®  and  the  simulation  began 
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ARRAY 


Fig.  D.4.  Variation  of  two-emitter  geometry  to  simulate 
dependence  for  thinned  linear  array . 


LOS 


^th  the  unfixed  signal  0.32  bea.«idths  ot  7.33“  away  from  the  anhiguous 
angle.  The  separation  was  reduted  in  3  dB  steps  heatwidth  separation  an 
then  to  a  aero  beamwldth  separation  (exact  coincidence  with  the  ambiguous 
angle).  The  second  input  variable  was  the  signal-to-nolse  ratio, 
consisted  of  100  "looks"  at  the  array  output,  that  Is,  100  data  vectors  whic 
were  processed  to  form  a  single  covariance  matrix.  After  the  ang  e 
estimation  algorithm  produced  a  set  of  direction  estimates  for  the  signals, 
the  power  estimation  algorithm  generated  a  power  estimate  for  direction.  For 
each  slgnal-to-nolse  ratio  (0  to  50  dB  in  10  dB  steps)  at  each  separation 
angle,  100  trials  were  conducted  in  order  to  compile  power  error  statistics 

for  each  solution  method. 

The  results  of  the  simulation  may  be  seen  in  Figs.  D.5  and  D.  . 
are  two  curves  for  each  separation,  one  for  each  signal,  and  the  area  between 
them  IS  shaded.  In  Fig.  0.5,  the  bottom  curve  (brown)  represents  the  best 
performance  of  the  direct  solution  method.  For  separations  greater  man  .32 
BB,  this  curve  approximately  represents  the  lower  error  bound.  When  the 
separation  decreases,  however,  representative  curves  break  away  from  the 
lower  bound  at  progressively  higher  signal-to-nolse  ratios.  Finally  the 

direct  method  breaks  down  completely  at  zero  separation. 

Now.  the  least-squares  performance  curves  In  Fig.  D. 6  may  be  compared 
with  the  direct  solution.  The  thick  curve  includes  all  of  the  separations  of 
Fig.  TI.5,  including  zero  separation.  The  curves  in  Fig.  D.6  are  all  almost 
exactly  coincident  with  the  lower  error  bound  for  the  direct  solution. 

regardless  of  the  separation  from  the  ambiguity. 

in  s«ry,  theory  and  experiments  show  that  the  direct  solution  power 
estimator  can  fail  for  non-uniform  arrays.  Whereas  the  traditional 
direction  finding  problem  has  been  resolving  closely-spaced  signals,  it  has 
been  shown  that  wide  angles  are  also  a  problem  when  non-uniform  arrays  are 
used.-  The  least-squares  power  estimator  is  able  to  identify  and  reject 
spurious  signals  In  situations  where  the  direct  solution  breaks  down  due  to 
ambiguities.  Experimental  data  shows  that  least-squares  performance  remains 
constant  independent  of  the  signal  separation  fro.  an  ambiguity.  Therefore, 
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RMS  POWER  ERROR 


Fig.  D.5. 

approach 


Comparison  of  direct  solution  power 
linear  dependence. 


estimation  error  as  emitters 


RMS  POWER  ERROR 


a  processing  system  using  this  power  estimation  technique  will  be  able  to 
overcome  the  ambiguity  problem  of  thinned  arrays  while  allowing  accuracy 
comparable  to  uniform  arrays. 
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Appendix  E 

SAMPLING  FROM  THE  WISHART  DISTRIBUTION 


1.  THE  COMPLEX  WISHART  DISTRIBUTION 

Let  r  be  an  M-dimensional,  complex  (circular)  Gaussian  vector  with  zero 
mean  and  covariance  R.  Given  K  observations  {r(k) |k=l, • • • .^}>  the  sample  co- 
variance  matrix 

K 

R  =  F  I  r(k)r®(k) 

^  k=l 


is  a  sufficient  statistic  for  inferring  parameters  of  R.  However,  when  K  is 
less  than  M,  the  sample  covariance  matrix  has  less  than  full  rank  and  is 
therefore  singular.  To  avoid  unnecessary  complications,  the  analysis  pre¬ 
sented  here  is  based  on  the  simplifying  assumption 

K  >  M 

Under  the  above  conditions,  the  probability  distribution  of  the  elements  of  a 
sample  covariance  matrix  is  completely  specified  by  the  complex  Wishart  dis¬ 
tribution  [31].  In  particular,  the  joint  probability  density  function  of  the 
elements  of 


A  =  K  R 


(E.l) 


may  be  written  as 


p(A) 


{  -Tr  (r”^A)  } 


The  normalization  constant  c(K,M)  is  given  by 


(E.2) 
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r(K)  ...  ra-M+i) 


c(K,M)  =  IT 


where  r(v)  is  the  gamma  function  [45] 


r( v)  =  /  e"^  dx  ;  V  >_  1 


Since  r(l)  =  1,  the  recursive  property  of  the  gamma  function 


r(vH)  =  V  r(v) 


(E.3) 


establishes  the  familiar  factorial  formula 

r(n)  =  (n-1)! 

for  any  positive  integer  n. 

The  probability  density  function  (pdf)  in  (E.2)  Is  defined  over  the  do¬ 
main  of  non-negative  definite  Hermltian  matrices.  While  theoretically  quite 
elegant,  this  formula  is  awkward  to  use  in  many  practical  calculations.  The 
fundamental  problem  is  that  the  elements  of  a  sample  covariance  matrix  are 
not  statistically  Independent.  Even  in  the  special  case  where  R  is  the  iden¬ 
tity  and  K=M,  the  restriction  of  (E.2)  to  a  semi -definite  domain  complicates 

matters  considerably. 

Fortunately,  the  Cholesky  decomposition  of  A  leads  to  a  much  more  useful 
representation.  Thus,  consider  the  class  of  (MxM)  upper  triangular  complex 
matrices  with  (real)  positive  elements  along  the  main  diagonal.  The  elements 
above  the  main  diagonal  are  unrestricted,  whereas  the  elements  below  the  main 
diagonal  are  identically  zero.  For  every  positive  definite  Hermltian  matrix 
A  defined  by  (E.l)  there  exists  a  unique  upper  triangular  matrix  U  such  that 


A  =  U^U 


112 


Goodman  131!  has  shown  that  the  joint  pdf  of  the  elements  of  U  Is  generally 
of  the  form 


However,  when  E  Is  the  identity  matrlm  I.  the  joint  pdf  of  the  elements  of  U 
may  be  empressed  as  the  product  of  the  marginal  densities  of  the  individual 
elements.  This  observation  stems  from  the  fact  that 


Tr  (U^J)  =  I 
m,n 

The  implication  is  that  in  the  special  case  E-I  the  elements  of  D  are  statis¬ 
tically  independent!  In  order  to  establish  this  Important  result,  we  first 

recognize  that 

(z)  =  (I/tt)  exp  {-|z|^}  ;  complex  z,  m  <  n  (E.5) 

mn 

is  the  pdf  of  a  (complex)  circular  Gaussian  random  variable  with  zero  mean 
and  unit  variance.  Setting  R=I  in  (E.4)  and  integrating  out  the  M(M-l)/2 
complex  Gaussian  terms  leaves  the  pdf  of  the  main  diagonal  of  U,  i.e., 

^  2  -2K-(2m-l)  \ 

p(TJli,  •••>  ~  ^  r(K-T3H-f)  nnn 

m— 1 


Evidently,  the  mth  diagonal  element  has  the  marginal  density 
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(E.6) 


Pit 

Tfirn 


2  „2(K-in+l)-l 

r(K-iiri-l)’' 


exp 


{- } 


,  X  >  0 


This  pdf  is  closely  related  to  the  chi-square  distribution  with  2(K-in+l)  de 
grees  of  freedom.  Thus,  consider  the  random  variable 


^n  = 


2 

^2n 


(E.7) 


obtained  by  simply  scaling  a  chi-square  random  variable  with  2n  degrees  of 
freedom.  It  can  be  shown  that  the  pdf  of  (E.7)  is  a  gamma  density,  i.e., 

Py  x""^  exp  {  -X  }  ,  X  >  0 


The  general  form  of  the  gamma  density  [46]  is  obtained  by  replacing  n  with  a 
real  number  1  and  introducing  a  positive  scale  factor  a,  i.e., 

p  (x)  =  a  p  (ax) 

V 


=  x'^^  exp  {  -ax  }  ,  X  >  0  .  (E.8) 

The  gamma  variate  (i.e.,  random  variable)  obtained  by  setting  a=l  in  (E.8)  is 
referred  to  here  as  the  vth  order  gamma  variate.  It  follows  from  the  recur¬ 
sive  property  of  the  gamma  function  (E.3)  that  the  expected  value  of  the  vth 

order  gamma  variate  is  v. 

The  pdf  of  the  square  root  of  the  vth  order  gamma  variate  is 
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p  (X)  -  2  ,  p  (x^) 

Y  V 

V 


2  2v-l 

"  r(v) 


exp 


,  X  >  0 


Comparing  this  pdf  with  (E.6),  we  conclude  that  the  diagonal  elements  of  U 
are  square  roots  of  integer  order  gamma  variates.  Specifically,  we  have 


^mm 

2.  APPLICATION  TO  MONTE  CARLO  SIMOLATIONS 

Sample  covariance  matrices  can  be  generated  directly  from  the  Wishart 
distribution  as  follows.  Given  an  arbitrary  covariance  matrix  R,  consider 
any  convenient  "square-root"  decomposition 

R  =  LL®  . 

The  observed  vector  r  may  be  interpreted  as  the  linear  transformation 


r  =  Lw 

of  a  normalized  vector  w  with  covariance 
E  {ww^ }  ~  ^ 

It  follows  that  any  sample  (covariance  matrix)  of  r  is  statistically  equiva¬ 
lent  to  a  linear  transformation  of  a  corresponding  sample  (covariance  matrix) 
of  w,  i.e.. 


R  =  LWL 


H 
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Since 


E  {w}  =  I 

the  argument  presented  in  the  previous  section  applies,  and  hence  we  may  con- 
struct 

K  W  = 

from  the  upper  triangular  matrix  U  with  statistically  independent  elements 
described  by  (E.5)  and  (E.6).  If  L  is  the  (lower)  Cholesky  factor  of  R,  then 

K  R  = 

=  (LU®)  (LD®)^ 

is  obtained  as  the  product  of  its  Cholesky  factors. 

The  complex  Gaussian  elements  of  U  (above  the  main  diagonal)  can  be  gen¬ 
erated  using  any  one  of  a  number  of  well-known  techniques.  Perhaps  the  most 
elegant  approach  is  to  compute 

z  =  (-In  u)^^^  exp  {  i2w  } 

where  u  and  v  are  statistically  Independent  random  variables  uniformly  dis¬ 
tributed  over  the  semi-open  Interval  (0,1].  This  method  is  based  on  the  fact 
that  -In  u  is  an  exponential  random  variable  with  unit  mean  (i.e.,  a  first 
order  gamma  variate).  Moreover,  nth  order  gamma  variates  can  be  generated  by 
summing  n  jointly  independent  exponential  variates.  Thus,  the  remaining  ele¬ 
ments  of  U  (on  the  main  diagonal)  could  be  obtained  from  expressions  of  the 

form 
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K-nH-1  1  /o 

X=  I  {-lnu(k)]'^^ 

where  the  {u(k) }  are  uniformly  distributed  over  (0,1]  and  jointly 
Independent.  However,  a  Fortran  subroutine  is  widely  available  [47]  that 
calculates,  gamma  variates  much  more  efficiently,  particularly  for  large 
values  of  K.  This  algorithm  is  based  on  the  method  of  acceptance/rejection 

testing  [48]. 
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APPENDIX  F 

DIRECTION  FINDING  EXPERHIENTS 


The  results  of  the  Monte  Carlo  experiments  described  in  Section  IV. A  are 
Included  in  their  entirety.  Four  different  pereaeter  yarletious  are 

explored : 

(1)  Signal-to-Interference  Ratio  Group 
A.  0  dB  B.  -10  dB 

(2)  Emitter  Separation  Group 

A.  0.1  Beamwidths  B.  0.2  Beamwidths  C.  0.4  Beamwidths 

(3)  Number  of  Looks  Group 
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10  B.  100 

C.  1000 
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Spectral  B. 

Root 
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ROOT  AAR 

MEM 

ROOT  MEM 
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ROOT  MLM 

MUSIC 

ROOT  MUSIC 

TNA 

Using  these  group  designations,  the  36  summary  plots  for  the  direction 
.  finding  experiments  can  be  indexed  as  follows : 
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Appendix  G 

ADAPTIVE  LISTENING  EXPERIMENTS 


The  results  of  the  Monte  Carlo  experiments  described  in  Section  IV. B. 
are  included  here  in  their  entirety.  The  curves  in  these  figures  all 
represent  loci  of  10  dB  output  signal-to-interference-plus-noise  ratio  for  a 
two-emitter  problem  with  the  desired  signal  10  dB  weaker  than  the 
interferer.  The  results  are  divided  into  three  sections:  (1)  Model 
Covariance  Method  Experiments  (pp.  G-3  -  G-14);  (2)  Projection  Method 

Experiments  (pp.  G-15  -  G-26);  and  (3)  Ef  f ects  .  of  Likelihood  Ratio  Test  for 
TiDSIC  (pp.  G-27  -  G-30). 

In  the  first  two  sections,  three  different  parameter  variations  are 
explored: 


(1)  DF  Method  Group 


A.  ARM*  B. 

ROOTMEM 

AAR 

ROOTTNA 

MEM 

ROOTMLM* 

TNA 

ROOTMDSIC^ 

MUSIC* 

MLM 

Snapshot  Group 

A.  20  B. 

100 

Calibration  Error  Group 

A.  0  dB  Amplitude  B. 

0.05  dB 

0®  Phase 

0.5® 

g  these  group  designations. 

the  24 

C. 


0.5  dB 
5.0® 


Covariance  Method  and  Projection  Nulling  Method  can  be  indexed  as  follows: 


*Note  that  in  the  projection  method  evaluation  only  the  algorithmns 
Indicated  were  compared,  since  the  results  tended  to  be  identical  with  those 
obtained  for  the  model  covariance  method. 
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The  remaining  four  plots  (G-25  to  -28)  illustrate  the  impact  of  the 
alternative  likelihood  ratio  tests  for  driving  MUSIC,  and  Root  MUSIC,  as 
described  in  Section  III  .A  and  Appendix  C.  It  can  be  seen  from  these  results 
that,  even  for  only  20  array  snapshots,  the  choice  of  LRT  is  not  important 
for  Signals  separated  by  less  the  0.4  beamvidths  with  spectral  MUSIC  or  less 

than  0.2  beamwldths  for  Root  MUSIC.  Moreover,  even  in  the  worst  case,  the 

Impact  of  the  choice  of  LRT  is  only  seen  to  be  at  most  a  couple  of  dB  in 

required  array  SNR.  Thus,  we  conclude  that  the  choice  of  LRT  is  more 

important  as  a  factor  in  direction-finding  performance  than  it  is  for 
adaptive  list eni ng . 
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